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Chapter  1 

Introduction 

Pressure  waves,  detected  by  an  array  of  receivers,  can  be  analyzed  to 
determine  the  location  of  the  acoustic  source,  or  the  location  of  objects  which 
the  waves  encountered  along  their  path.  This  thesis  studies  those  methods  of 
analysis  which  can  provide  particularly  precise  estimates  of  the  locations.  Such 
methods  are  typically  referred  to  as  “super-directivity,”  “super-resolution,”  or 
“optimum”  processors. 

Subsequent  sections  of  this  introduction  will  describe  the  questions  which 
motivated  the  study,  list  the  assumptions  which  will  be  used  in  the  analysis  of 
the  pressure  waves,  state  the  scope  of  the  study,  and  then  conclude  with  a  spec¬ 
ification  of  some  of  the  symbols  and  relationships  used  in  subsequent  chapters. 

1.1  Motivation  for  the  Study 

Theory  and  experience  in  the  processing  of  the  signals  obtained  from  an 
array  of  receivers  has  accumulated  steadily,  with  one  of  the  earliest  investigations 
being  reported  by  Oseen1  in  1922.  Because  arrays  and  pressure  waves  constitute 
a  very  general  phenomenon,  related  discoveries  have  been  made  in  many  different 
disciplines.  This  has  resulted,  as  pointed  out  by  Miller,2  in  the  same  fundamental 
principle  sometimes  being  discovered  more  than  once. 

As  the  basic  concepts  of  existing  methods  were  studied,  it  became  clear 
that,  because  of  the  inter-disciplinary  nature  of  the  problem,  it  was  possible  that 
discoveries  in  one  field  might  have  been  overlooked  by  workers  in  another.  Thus, 
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one  motivation  for  this  study  was  to  consider  combinations  of  methods,  from 
perhaps  different  disciplines,  which  could  be  combined  in  new  and  productive 
ways.  In  particular,  the  approach  of  using  a  ^fictitious  source,”  discussed  in 
Chapter  2,  is  combined  with  the  concept  of  a  matched  filter  in  Chapter  4  to 
develop  two  new  imaging  algorithms  called  the  pattern-match  and  mismatch 
methods. 

Most  previous  work  has  concentrated  quite  naturally  on  the  pragmatic 
cases  where  one  wishes  to  locate  a  sound  source  in  noise.  Some  of  the  more  imagi¬ 
native  methods  are  non-linear  or  iterative  and  so  defy  both  traditional  closed-form 
analysis  and  a  straightforward  conceptual  understanding.  The  author’s  interest 
was  thus  piqued  by  the  following  question:  which  methods  are  primarily  for  noise 
suppression,  and  which  are  utilizing  some  new  insight  into  the  underlying  wave 
phenomena?  In  addition,  we  wished  to  know  if  these  methods,  which  were  de¬ 
veloped  for  waves  from  sound  sources,  could  also  be  applied  to  waves  from  sound 
scatterers.  The  answer  to  these  questions  could  be  useful  in  matching  the  various 
methods  to  a  given  problem. 

Finally,  several  methods  touched  upon  the  idea  of  using  a  swept-frequency 
transmitted  signal.  Indeed,  linear  FM  sweeps  have  long  been  used.  However, 
many  methods  were  built  upon  the  assumption  that  the  object  which  was  en¬ 
countered  by  the  pressure  wave  was  a  point-reflector;  such  an  assumption,  it 
seemed,  may  cause  potentially  useful  information  to  be  lost.  A  final  motivation, 
then,  was  to  explore  a  new  frequency-swept  method  which  was  built  upon  the 
assumption  of  frequency-dependent  scattering  from  these  objects. 
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1.2  Assumptions  About  the  Propagation  Path 

The  modeling  of  the  complex  phenomenon  of  a  realistic  environment  is 
certainly  a  worthy  and  demanding  pursuit;  however,  in  this  study,  a  far  simpler 
situation  will  be  modeled  in  order  that  the  fundamental  principles  can  be  better 
understood  and  simulated. 

This  study  will  generally  assume  that  the  wave  phenomenon  occurs  in  an 
isotropic,  unbounded  medium  wherein  the  speed  of  sound  is  known  and  damping 
is  negligible.  Wave  propagation  will  be  modeled  as  spherical  or  planar  compres- 
sional  waves  obeying  the  linear  wave  equation. 

In  most  cases,  the  source  of  the  pressure  waves  will  be  an  omni-directional 
transmitter  whose  location  is  known  and  which  can  broadcast  acoustic  energy 
of  any  desired  temporal  pattern.  This  pattern  will  typically  be  a  windowed  sine 
wave  at  a  given  frequency,  or  a  series  of  such  sine  waves  which  step  through  a 
set  of  frequencies.  Amplitudes  will  be  arbitrary  but  shall  be  low  enough  that  the 
assumption  of  linear  propagation  is  valid.  These  assumptions  effectively  eliminate 
any  question  about  the  location  of  the  source  of  the  pressure  waves;  thus,  this 
study  will  focus  on  discovering  the  precise  location  of  objects  which  the  waves 
encountered  along  their  path.  However,  many  methods,  which  were  originally 
developed  to  find  the  source  of  acoustic  energy,  can  be  recast  to  find  these  objects. 

The  detection  of  the  pressure  waves  will  be  simulated  as  an  equidistant 
line  array  of  omni-directional,  linear  receivers.  Except  for  certain  non-linear 
spectral  estimators,  the  relationship  between  a  one-  and  two-dimensional  array 
is  straightforward  and  so  this  study  will  be  limited  to  the  one-dimensional  array. 
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It  will  be  assumed  that  the  objects  which  are  encountered  by  the  waves 
have  a  velocity  of  zero.  When  more  than  one  object  is  encountered,  it  will  be 
assumed  that  the  Born  approximation  for  weak  scattering  can  be  applied.  Most 
objects  will  be  modeled  as  point  scatterers,  spheres,  or  disks. 

1.3  Scope  of  the  Thesis 

We  wish  to  find  methods,  given  the  assumptions  described  in  Section  1.2, 
for  analyzing  samples  Rj  of  the  amplitude  of  an  pressure  wave  which  yield,  with 
the  highest  possible  precision,  the  location  of  the  object  or  objects  encountered 
by  the  pressure  wave.  The  samples  will  be  obtained  at  locations  X,  which  are 
equally  spaced  along  a  line.  In  the  pursuit  of  these  methods,  all  the  existing  basic 
and  high  resolution  methods  will  be  considered;  the  most  promising  concepts, 
along  with  some  new  approaches,  will  be  combined  to  form  several  new  acoustic 
imaging  methods  which  will  then  be  explored  through  simulations. 

1.4  Definition  of  the  Major  Terms 

The  path  of  the  pressure  wave  begins  at  the  transmitter  which  is  located 
at  position  [XE,ZE^J  as  shown  in  Figure  1.1.  The  pattern  for  a  windowed  sine 
wave  of  angular  frequency  ui  (=  2nf)  is  given  by 

(  sinful ),  if  0  <  t  <  T; 

m  =  \  .  (i.i) 

f  0,  otherwise 

where  the  period  or  duration  of  the  signal  is  T  and  the  time  t  is  continuous. 

The  wave  next  encounters  an  object  at  some  unknown  position  (X5,  Zs^j 
with  a  scattering  function  S  (<j>R\<j>E'\ .  For  a  point  scatterer,  5  is  a  constant 
(usually  1).  Otherwise,  it  is  a  function  of  two  angles  which  are  defined  relative  to 


scattering  object 


the  scatterer.  Given  that  the  pressure  wave  is  arriving  at  angle  4>E ,  the  scattering 
function  then  gives  the  amplitude  of  the  wave  which  is  reflected  into  the  angle 
<j>R.  The  scattering  function  will  be  simulated  in  most  cases  as  causing  no  phase 
shift,  i.e.,  it  is  a  real  rather  than  complex  function.  When  we  need  to  consider 
more  than  one  scattering  object,  each  object  will  be  at  position  with 

scattering  function  Sm  (tmltm)  1  <  m  <  NS  where  there  are  Ns  objects. 

The  pressure  wave  is  sampled  at  the  receivers  located  at  positions  (XR,  o) 
for  1  <  i  <  Nr  where  there  are  NR  receivers.  The  difference  between  any  two 
receiver  positions  (the  spatial  sampling  interval)  is  AX.  The  Z  coordinate  is 
always  zero  since  we  have  chosen  to  place  the  line  array  along  the  X-axis.  The 
amplitude  of  the  pressure  wave  )  are  sampled  at  times  f/  for  1  <  /  <  NT 

where  there  are  NT  time  samples.  The  difference  between  any  two  time  samples 
(the  temporal  sampling  interval)  is  At. 


The  signal  detected  at  each  receiver  can  be  represented  as  the  linear  sum 
of  the  pressure  reflected  from  the  one  or  more  objects: 

NS  i 

RiM  =  £  n - rT~  5™  (*£.  **)  ■ E  ('' "  T™)  '  (L2> 

-  =  'Jixf-xg)2  +(Z«)2;  (1.3) 

DTm  =  \/(XT  -  X*f  +  (ZT  -  Zif;  and  (1.4) 

Tfm  =  D,m  +  °2"  .  (1.5) 


The  time  of  flight  T[m  from  the  transmitter  to  the  ith  receiver  via  a  reflection 
from  the  mth  object  is  calculated  using  the  length  of  the  flight  path,  composed 
of  the  path  from  the  transmitter  to  an  object  (distance  DTm )  plus  the  path  from 
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the  object  to  a  receiver  (distance  Z),m),  for  waves  traveling  at  the  speed  of  sound 
c  where  c  =  A/. 

When  the  locations  of  the  objects  are  known,  they  can  be  combined  with 
the  location  of  the  transmitter  and  receivers  to  compute  the  scattering  angles  as 


±E  (XT-X ^ 

K  =  arctan 

(yS  _  yR ' 

my—  ~ 


,  and 


(1.6) 

(1.7) 


A  matrix  notation  will  be  used  in  some  cases  so  that  some  of  the  above 
quantities  can  be  expressed  more  compactly  as 

R(t{)  =  [Ri(tt)  R2(ti)  ...  and  (1.8) 

XR=[X?  X?  ...  X§]T  (1.9) 

where  T  denotes  the  transpose  operation  which  make  these  representations  col¬ 
umn  vectors.  When  there  is  only  one  object,  Eq.  (1.8)  becomes 


R(U) 


... 

...-5-L_S^)-E(t,-TZj\T 


(1.10) 


where  the  implicit  subscript  on  S  and  the  implicit  second  subscript  on  TF  is  1 
since  there  is  only  one  object. 
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Chapter  2 

Basic  Issues  in  Array  Processing 

Subsequent  chapters  will  investigate  how  various  processing  methods  can 
increase  a  receiver  array’s  resolution.  This  chapter  will  lay  the  groundwork  for 
those  investigations  by  examining  the  basic  approaches  and  issues. 

The  problem  being  addressed  by  this  study  has  been  approached  from 
several  points  of  view  in  various  fields  which  seem,  at  first  glance,  to  be  only 
vaguely  related.  This  does  not  seem  reasonable,  however,  since  the  underlying 
physics  is  the  same  in  each  case.  Therefore,  we  wish  to  roughly  categorize  these 
approaches,  and  then  seek  a  unified  approach  to  which  all  the  others  are  special 
cases.  We  expect  the  various  cases  to  differ  only  in  the  nature  of  their  assumptions 
or  their  requirement  of  a  priori  information. 

2.1  Holographic  Reconstruction 

In  the  traditional  study  of  the  propagation  of  pressure  waves,  the  inter¬ 
action  of  the  waves  with  an  object  is  categorized  as  diffraction.  The  problem  we 
are  addressing  in  this  study  is  then  categorized  as  inverse  diffraction.  In  particu¬ 
lar,  the  holographic  reconstruction  methods  use  inverse  diffraction  to  recreate  or 
reconstruct  the  pressure  field  in  the  vicinity  of  objects  which  are  either  scatterers 
or  are  themselves  acoustic  sources. 
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2.1.1  Holographic  Assumptions 

In  the  typical  three-dimensional  case,  the  receivers  are  equally  spaced 
across  a  planar  grid,  and  the  objects  are  assumed  to  be  located  on  a  plane  which 
is  parallel  to  the  receiver  plane.  For  this  reason,  most  holographic  studies  have 
concentrated  on  reconstructing  this  object  plane  which  is  then  interpreted  as  an 
“image”  of  the  objects.  The  image  can  be  formed  to  represent  the  amplitude  or 
phase  of  the  pressure,  the  vector  velocity,  the  vector  intensity,  or  other  quantities. 

The  interpretation  of  the  image  is  by  no  means  deterministic  -  it  depends 
upon  the  expected  nature  of  the  objects.  If  point  objects  are  expected,  each  point 
of  the  image  could  be  interpreted  as  the  amplitude  of  the  object  located  at  that 
point  where  amplitudes  which  fall  below  some  threshold  would  indicate  that  no 
object  is  present.  On  the  other  hand,  if  extended  objects  are  expected,  each  point 
of  the  image  could  be  interpreted  as  the  amplitude  of  the  pressure  scattered  or 
generated  at  that  location  on  the  object. 

The  earliest  holographic  work,  such  as  that  by  Graham3  and  Watson,4,5 
used  optical  methods  to  process  the  information.  This  background,  and  the 
limited  speed  of  the  computers  of  the  time,  seems  to  have  motivated  the  early 
users  of  digital  processing  to  make  the  Fraunhofer  or  Fresnel  approximation  as 
it  eliminated  the  evaluation  of  a  square  root  in  the  algorithm.  With  the  work 
of  IBM’s  VanRooy,6,7  digital  processing  eliminated  this  approximation  and  the 
approach  became,  therefore,  valid  for  both  the  near  and  far  fields.  VanRooy  was 
also  one  of  the  first  to  apply  the  Fast  Fourier  Transform  (FFT)  which  significantly 
reduced  the  computer  time  needed;  his  basic  algorithm  continues  to  be  used 
over  a  decade  later  for  holographic  reconstructions  in  studies  such  as  those  by 
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Stepanishen  and  Benjamin,8  Eschenberg  and  Hayek,9,10,11  and  Williams  et  a/..12 
2.1.2  The  Holographic  Equation 


A  variety  of  methods  can  be  used  to  develop  the  holographic  equation 
which  expresses  the  pressure  field  in  the  object  plane  as  a  function  of  the  measured 
field  in  the  receiver  or  hologram  plane.  Cohen13  used  the  Rayleigh  integral, 
expressed  as  a  convolution  in  the  space  domain  and  as  a  multiplication  in  the 
spatial  frequency  domain,  to  develop  his  method  of  processing.  Another  approach 
based  upon  a  Fourier  transform  of  the  wave  equation  is  shown  in  Appendix  A. 


One  of  the  earliest,  and  still  most  elegant,  developments  was  provided  by 
Goodman.14  In  the  two-dimensional  system  of  this  study,  the  Helmholtz  wave 
equation  for  harmonic  time  dependence  becomes 


V2R(x,  z)  +  k2R(x,  z )  =  +  k2R(x,  z)  =  0.  (2.1) 


The  signal  R(x,  z)  can  be  written  in  terms  of  an  inverse  Fourier  transform  as 


R(x,  z)  =  T  1  ( R(kx,z )) 

1  +°° 

=  — ;=  J  R{kx,  z)  ■  exp  {-j  kx  x)  dkx 

— OO 


(2.2) 


where  T  1  represents  the  inverse  Fourier  transform  while  T  and  R  represents 
the  forward  Fourier  transform  as  given  by 

T(R{x,  z)J  =  R(kx,  z)  =  — p=  J  R(x,  z)  ■  exp  (j  kx  x)  dx.  (2.3) 

^  — oo 


The  letter  j  denotes  the  square  root  of  minus  one.  These  definitions  of  the  forward 
and  inverse  seem  to  agree  with  the  majority  of  scientific  users  of  the  Fourier 
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transform13  but  are  reversed  from  some,  including  the  ones  used  in  some  previous 
work  by  the  author.11  However,  simulations  have  proven  that  these  definitions 
can  be  reversed  with  no  effect  on,  at  least,  the  reconstructed  pressure  amplitude. 
Upon  substitution  of  Eq.  (2.2)  into  Eq.  (2.1),  we  note  that  the  integrand  of  the 
first  term  becomes  — k 3  R(kx,z).  Using  the  relationship  between  the  components 
of  the  wavenumber, 

k'  =  k“  +  k~,  (2.4) 


we  can  combine  the  integrands  to  yield 


+oo  r 

!_  [ 


d2R(kx,z) 
dz 2 


+  k2  R(kx,z) 


exp  (—  j  kx  x)  dkx  —  0. 


(2.5) 


The  integral  can  only  be  zero  for  all  values  of  the  independent  variable  x  if  its 
integrand  is  zero;  thus,  we  can  write 


d2R{kx,z) 
dz 2 


+  k2  R(kx,z)  =  0 


(2.6) 


which  is  the  transformed  version  of  the  wave  equation.  One  solution  to  this 
differential  equation  is 


R{kx,z)  =  R(kx,0)  ■  exp(j  kz  z)  (2.7) 

where  R(kx,0)  is  often  written  as  R^ (kx)  to  emphasize  that  it  is  measured  in 
the  receiver  or  hologram  plane.  Finally,  we  can  substitute  Eq.  (2.7)  into  Eq.  (2.2) 
so  that  the  field  for  a  plane  at  any  location  a  can  be  written  as 

R{x,  z)  =  -j=  J  RH{kx)  ■  exp  (j  kz  z  -  j  kx  x)  dkx.  (2.8) 

—  OO 

It  is  sometimes  useful  to  refer  to  that  part  of  this  equation  which  represents 
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the  effects  of  propagation  through  space  as  the  Green’s  function  in  the  spatial 
frequency  domain: 

G(kz)  =  exp  (j  kz  z) .  (2.9) 

2.1.3  Discrete  Holography 

The  conversion  of  the  continuous  version  of  the  holography  equation, 
Eq.  (2.S),  to  its  discrete,  sampled  form  has  been  described  previously11,13  but 
is  repeated  here  for  completeness,  using  the  terminology  of  this  study.  Beginning 
with  the  forward  Fourier  transform  of  the  hologram,  Eq.  (2.3),  we  replace  the 
continuous  variables  with  the  discrete  variables: 

dx  — *•  AX , 

x  —>  i  AX, 

dkx  -*•  AKX  =  (2.10) 

L/ 

L  =  Nr  AX, 
kx  *  (j  AXx 

It  is  important  to  note  that  the  receiver  aperture  size  L  is  not  merely  the  dis¬ 
tance  between  the  receivers  at  the  ends  of  the  array  but  is  AX  larger  than  that 
distance.  The  signal  in  the  receiver  aperture,  RH(x),  is  now  sampled  at  the 
discrete  receiver  locations  XR  =  i  AX  and  so  could  be  written  RH(i  AX)  or 
just  Rf .  Likewise,  the  signal  in  the  spatial  frequency  domain  is  sampled  and  so 
RH (kx)  becomes  RI{ (q  A I(x)  or  just  Rr .  Substituting  into  Eq.  (2.3),  the  forward 
transform  becomes 

av  Nr 

=  -s=  E  R>  ■  “P  U  qAKxi&X). 

V2ir  ,=1 


(2.11) 
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With  similar  substitutions,  the  holographic  reconstruction,  Eq.  (2.S),  becomes 


a  K 


R{x,  z)  =  —j==  Y1  ~Rq  ■  exP  O'  kz  z  -  j  x  q  A Kx) 
v2tt  ,=1 


(2.12) 


ks  =  \A2  -  {q  A Kxy 


(2.13) 


Eq.  (2.12)  is  valid  for  any  value  of  FV s  argument,  the  continuous  variables  x  and 
2.  However,  if  this  is  to  be  implemented  as  an  FFT,  we  must  also  convert  x  to 
its  discrete  form  so  that  we  have 

A  K  nR 

R(i  AX ,  z)  =  Ri(z )  =  —J=  Y]  M  ■  exp  (j  kz  z  -  j  q  A Kx  i  AX) .  (2.14) 

?=:1 


In  many  studies,  the  indices  used  in  discrete  equations  are  arranged  so 
that  they,  as  well  as  the  discrete  values  of  x  and  fcz,  are  more  or  less  symmetrical 


about  zero.  That  is,  the  indices  are  transformed  so  that 


1  <i<  NJ 


1  <q<N* 


-\NR  <  i  <  \NR  -  1, 
-%NR  <q<  %NR-  1. 


(2.15) 


Note  that  since  0  is  included  in  the  shifted  domain,  we  must  subtract  1  from 
one  of  the  limits  so  that  the  total  number  of  points  continues  to  be  NR.  As  is 
obvious  from  Eq.  (2.11),  these  shifts  are  equivalent  to  adding  a  constant  value  to 
the  phase  of  the  results;  therefore,  we  can  choose  whichever  domain  we  prefer  as 
long  as  we  are  consistent. 


There  are  two  reasons  for  shifting  the  domain  of  q.  First,  a  symmetrical 
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domain  reflects  the  fact  that  both  positive  and  negative  spatial  frequencies  are 
involved.  That  is,  index  q  refers  to  positive  and  negative  spatial  frequencies  as 
follows: 


f  1  <  q  <  \NR  +  1, 

unshifted:  < 

{  ±NR  +  2  <q<NR , 


and 


shifted: 


(  0  <  q  <  \XR, 

+  l  <  q  <  — 1, 


positive  spatial  frequencies; 
negative  spatial  frequencies; 


positive  spatial  frequencies; 
negative  spatial  frequencies. 


(2.16) 


(2.17) 


Second,  most  computer  implementations  of  the  FFT  are  arranged  so  that  the 
spatial  frequencies  are  returned  in  the  order  of  the  shifted  domain.  That  is, 
the  first  element  of  the  returned  data  array  corresponds  to  q  =  —\NR.  This, 
however,  requires  only  that  we  be  careful  when  combining  analytically  computed 
terms  with  those  computed  by  an  FFT. 


There  does  not  seem  to  be  any  particular  reason  for  shifting  the  domain 
of  x  except  for  perhaps  tradition  and  to  mimic  the  shifting  applied  to  q. 


2.1.4  Definition  of  Phase 


The  sign  of  the  exponent  in  Eq.  (2.7)  implies  a  choice  which  should  be 
noted  and  maintained,  for  consistency,  in  other  equations  in  this  study.  The 
general  solution  to  the  wave  equation  for  plane  waves  is  of  the  form 

<f>  =  A  ejwt~]kT  +  B  e]wt+}kT  (2.18) 


where  A  and  B  are  complex  numbers  and  r  is  in  the  direction  of  propagation.  As 
noted  by  Skudrzyk,15  the  choice  of  sign  for  jut  is  arbitrary;  however,  by  using  jut 
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instead  of  —  jut,  a  similarity  with  electrical  equations  is  maintained.  Though  this 
representation  is  for  plane  waves,  the  following  discussion  of  phase  and  direction 
will  apply  to  spherical,  cylindrical,  or  any  other  representation. 

Both  terms  will  be  needed  to  express  steady-state  conditions  or  situations 
involving  reflections.  However,  simple  propagation  can  be  simplified  by  using 
only  one  of  the  terms.  The  choice  of  term  can  be  made  by  observing  that  the 
first  corresponds  to  wavefronts  moving  so  that  r  is  increasing  while  the  second 
corresponds  to  wavefronts  where  r  is  decreasing.  In  most  cases  we  are  interested 
in  waves  which  are  traveling  away  from  their  source,  where  r  is  increasing,  and 
so  quite  naturally  choose  to  use  only  the  first  term.  However,  where  we  are 
interested  in  waves  traveling  back  towards  a  source  (such  as  Eq.  (2.7)),  we  will 
need  to  use  the  second  term. 

It  is  important  to  note  that,  because  we  have  chosen  the  first  term  of 
Eq.  (2.18),  the  phase  is  given  by  —hr  instead  of  kr.  Pressure  waves  arriving  at 
a  linear  array  will  take  longer  to  arrive  at  the  distant  elements  than  to  arrive  at 
the  near  elements.  While  the  absolute  value  of  the  pha=*>  will  be  greater  at  the 
distant  elements  (because  of  the  r),  the  actual  numerical  value  will  be  smaller 
(because  of  the  minus  sign).  Thus,  the  simulations  will  model  the  phase  received 
at  the  distant  elements  as  being  less  than  the  phase  received  at  the  near  elements. 
Although  this  definition  of  phase  is  common  among  theoretical  studies,  it  is  the 
opposite  from  that  derived  from  a  simple  experimental  system  where  the  greater 
time  delay  experienced  by  the  distant  elements  of  the  array  would  have  been 
recorded  as  a  greater  phase  delay. 
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2.1.5  Evanescent  Waves 

As  pointed  out  by  Goodman,1'1  each  Fourier  component  of  the  image  can 
be  viewed  as  a  plane  wave  propagating  with  a  direction  cosine  given  by 


a  =  A  kx . 


(2.19) 


These  waves  have  an  imaginary  direction  cosine  and  are  exponentially  attenuated 
as  they  leave  their  source.  As  a  result,  when  they  are  reconstructed  by  Eq.  (2.21), 
they  should  be  exponentially  amplified.  The  square  roots  in  Eq.  (2.20)  each  have 
a  positive  and  a  negative  solution;  we  have  chosen  the  negative  root  for  the 
evanescent  waves  so  that  the  mathematics  properly  produce  this  exponential 
amplification. 


Cohen13  calls  the  above  algorithm  the  backward  tracer  since  it  traces  all 
plane  waves  back  to  the  source.  In  most  experimental  systems  where  the  hologram 
is  recorded  several  wavelengths  or  more  from  the  objects,  the  value  of  Rh{^x) 
for  evanescent  waves  will  represent  only  noise,  numerical  roundoff,  and  aperture 
effects  since  these  waves  will  have  decayed  to  almost  nothing  after  several  wave- 
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lengths.  The  backward  tracer  will  therefore  exponentially  amplify  only  noise, 
producing  either  poor  or  utterly  useless  images.  As  an  alternative,  Cohen  used 
what  he  called  the  backward  propagator  formed  by  defining 


kz 


k  x  ^  k ; 
kz  >  k. 


(2.22) 


By  choosing  the  positive  square  root  for  the  kx  >  k  case,  the  evanescent  waves 
are  reconstructed  using  an  exponential  with  a  negative,  rather  than  positive,  real 
power;  that  is,  the  evanescent  waves  are  filtered  out.  This  approach  corresponds 
to  treating  each  point  on  the  hologram  as  being  itself  a  source  with  an  amplitude 
and  phase  given  by  R[i(x );  thus,  the  reconstructed  image  is  due  to  the  propaga¬ 
tion  of  these  artificial  sources  back  to  the  plane  of  the  actual  sources.  Since  the 
propagation  is  viewed  as  being  away  from  the  artificial  sources,  the  evanescent 
waves  are  exponentially  attenuated. 


When  working  in  the  extreme  nearfield,  evanescent  waves  can  be  recorded 
and  may  be  necessary  in  the  reconstruction  of  the  details  of  the  scattering  surface 
or  the  source  of  acoustic  energy.  Even  so,  the  signal-to-noise  ratio  for  these 
components  may  be  quite  low  and  thus  some  sort  of  filtering  is  appropriate.  For 
a  particular  experimental  holographic  system,11  it  was  noted  that  evanescent 
waves  which  had  decayed  by  12dB  or  more  could  not  be  reliably  recorded  and  so 
should  not  be  reconstructed.  This  approach  was  then  used  to  construct  a  filtered 
version  of  the  Green’s  function: 

(l  -  \  exp  (-6.666(1  -  kx/kc)j^,  kx  <  kc , 

kx  =  kc, 

2  exp (-6.666  (kx/kc  -  1)),  kx  >  kc , 


GHz) 

exp  (j  kz  z) 


(2.23) 


i 
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where  the  factor  of  — C.666  was  chosen  to  yield  a  relatively  sharp  slope,  and  the 
term  kc  is  chosen  so  that  G p  yields  at  most  a  12db  amplification. 

2.1.6  Some  Limitations  to  the  Holographic  Approach 

The  Helmholtz  wave  equation,  upon  which  the  holographic  equation  is 
built,  is  valid  only  for  isotropic  media  without  obstructions  or  sources.  If  there 
were,  say,  an  infinite,  perfectly  reflecting  plane  surface  in  the  vicinity  of  the  re¬ 
ceiver.  the  Green's  function,  Eq.  (2.9),  would  need  an  additional  term  to  represent 
this  reflection.  More  complicated  situations  with  multiple  reflectors  would  lead 
to  extremely  awkward  and  complicated  Green’s  functions.  It  may  become  impos¬ 
sible  to  write  the  Green’s  function  in  the  spatial  frequency  domain  as  is  required 
by  Eq.  (2.S).  As  we  shall  see  later,  other  approaches  can  provide  a  much  more 
straightforward  method  of  handling  these  situations. 

In  some  cases,  one  may  have  a  priori  knowledge  of  the  objects  reflecting  or 
generating  the  pressure  waves.  As  a  general  principle,  such  knowledge  should  be 
used  if  at  all  possible  to  enhance  the  imaging  process.  However,  there  is  no  direct 
way  to  do  this  using  the  holographic  approach  since  it  describes  only  pressure 
waves  in  free  space  rather  than  sources.  Technically  speaking,  the  holographic 
reconstruction  can  only  generate  valid  images  up  to,  but  not  including,  the  plane 
of  the  objects. 

One  of  the  great  advantages  of  the  holographic  approach,  as  formulated 
in  Eq.  (2.8),  is  that  it  consists  mostly  of  Fourier  transforms  which  can  be  im¬ 
plemented  efficiently  on  a  computer  using  the  FFT  algorithm.  However,  this 
efficiency  is  lost  and  the  equations  can  become  awkward  when  the  line  (or  plane) 
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of  the  receivers  and  the  line  (or  plane)  of  the  objects  are  not  parallel.  Approxi¬ 
mations  to  non-planar  objects  can  be  made  by  reconstructing  a  series  of  closely 
spaced  lines  (or  planes)  though  this  still  increases  the  amount  of  computations. 
In  some  cases,  a  space-domain  convolution,  as  investigated  by  Cohen,16  may  be 
more  efficient  even  though  it  does  not  use  the  FFT  at  all. 

The  image  also  becomes  ambiguous  when  there  are  several  discrete  objects 
which  do  not  lie  in  a  plane  but  are  instead  spread  out  over  two  (or  three)  dimen¬ 
sions.  There  is  no  closed-form  method  known  to  detect  or  reliably  analyze  this 
situation  although  Lang17  had  some  success  by  using  pulses  and  a  window  in  time 
that  accepted  only  signals  from  a  given  distance.  Narasimhan  et  a/.18  proposed 
two  methods  for  determining  the  range  for  the  special  case  of  a  plane  reflector 
illuminated  by  plane  waves.  Powers  and  Mueller19  developed  an  algorithm  that 
searched  for  planes  exhibiting  peaks  in  the  reconstructed  pressure  field;  planes 
with  the  highest  peaks  were  then  considered  to  contain  the  “true”  objects  while 
the  other  planes  were  considered  to  be  merely  out-of-focus  versions  of  the  true 
objects. 

A  similar  ambiguity  occurs  even  when  there  is  only  one  object  if  that 
object  itself  extends  through  several  planes.  For  example,  Lang20  showed  the 
reconstruction  of  two  planes  of  a  cone  whose  tip  is  pointed  towards  the  receiver 
array:  one  plane  at  the  tip,  and  one  plane  at  the  back  face  of  the  cone.  In  general, 
it  was  difficult  to  distinguish  between  the  phenomena  from  these  two  planes  — 
for  example,  the  circle  in  the  plane  of  the  back  face  could  have  either  been  the 
effects  of  the  backface  itself,  or  an  out-of-focus  image  of  the  tip. 


One  way  to  summarize  these  last  few  limitations  is  to  note  that,  since  it 
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is  basically  a  single-frequency  algorithm,  the  holographic  approach  has  no  range 
information  and  is  therefore  ill-suited  to  any  situation  where  range  information  is 
important.  Some  studies21, 22,23,24  have  formed  images  by  combining  the  images 
made  at  different  frequencies;  none,  however,  seem  to  have  really  overcome  the 
range-focusing  problem  within  the  realm  of  holography.  It  seems  that  other 
approaches,  such  as  will  be  discussed  later  in  this  chapter,  are  needed. 

2.2  Beamforming 

Perhaps  the  most  traditional  method  of  processing  the  output  of  an  array 
is  beamforming.  It  differs  from  holography  in  that  only  the  direction  (and  perhaps 
the  range)  of  the  objects  are  of  interest  —  it  is  not  used  to  form  an  image.  The 
central  concept  is  to  choose  weights  Ai  and  time  delays  T{  to  reinforce  pressure 
waves  arriving  from  a  given  direction  0: 

Nr 

B(e,tl)  =  Y,Al-R,(ti~Ti).  (2.24) 

i=i 

A  great  many  variations  can  be  built  upon  this  basic  formula. 

2.2.1  Bartlett  Beamformer 

One  of  the  most  widely  used  variations  is  the  so-called  Bartlett  beam- 
former.25  It  assumes  that  the  signals  consists  solely  of  plane  waves  of  a  given 
frequency  /,  and  that  differences  in  time  can  be  adequately  represented  by  phase. 
Since  differences  in  time  which  are  greater  than  one  period  (1//  seconds)  cannot 
be  unambiguously  represented  as  a  phase,  this  assumes  that  either  the  time  dif¬ 
ferences  are  sufficiently  small,  or  that  the  ambiguity  is  insignificant.  The  received 
signal  R,  can  then  be  expressed  as  a  complex  number  which  is  often  calculated 
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by  taking  the  Fourier  transform  of  the  received  signal’s  time  history.  The  time 
dependence  f/  in  the  beam  output  is  also  replaced  with  a  phase  representation, 
i.e.,  B  becomes  complex.  Finally,  we  can  express  the  time  delay  r,  as  the  phase 
of  the  now  complex  weights  to  write 


N 


R 


B{6 )  =  Y^A*-Rt  =  A'  R 


(2.25) 


i=i 


where 

A  =  [Ai  A'i  ...  Ax]T, 

A{  =  exp  kXf1  sin(0))  , 


(2.26) 


the  asterisk  (*)  denotes  the  complex  conjugate  operation,  the  T  denotes  the  trans¬ 
pose  operation,  and  the  dagger  (f)  denotes  both.  This  formulation  sets  the  point 
of  zero  phase,  which  is  arbitrary,  to  be  at  the  origin  of  the  coordinate  system. 
Often  the  power  in  the  beam  is  of  interest  and  is  computed  from  Eq.  (2.24)  as 

B2{0,t,)  =  BB* 


=  (Af  R)  ■  (Af  R)* 
=  A]  RR*  A 


(2.27) 


=  A^A 


where  3?  is  the  spatial  cross-correlation  matrix.  This  is  a  zero  time-delay  corre¬ 
lation  since  no  delays,  beyond  the  ones  naturally  occurring  due  to  the  time  of 
flight,  are  introduced.  For  example,  in  its  simplest  form,  3?;/,  =  Ri(ti)  ■  )  for 

any  two  receivers  i  and  h  including  i  =  h. 
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2.2.2  Limitations  of  the  Bartlett  Acoustical  Model 


To  review,  the  temporal  spectrum  can  be  thought  of  as  an  attempt  to 
replace  the  original  time  series  with  a  sum  of  sine  waves.  Thus,  this  spectrum 
has  difficulty  with  time  series  that  are  inherently  not  sine  waves.  For  example,  a 
step  function  cannot  be  represented;  even  if  the  spectrum  extends  to  infinity,  a 
small  overshoot,  called  Gibb’s  phenomenon,26  will  remain. 

In  an  analogous  manner,  the  Bartlett  beamformer  attempts  to  find  sine 
waves  across  an  equally  spaced  linear  array  where  each  sine  corresponds  to  a  plane 
wave  arriving  from  a  specific  direction.  Because  this  beamformer  is  modeling  the 
received  field  as  plane  waves,  it  will  have  difficulty  with  pressure  waves  which 
are  inherently  not  plane  waves,  such  as  spherical  waves.  This  problem  is  often 
referred  to  as  “model  mismatch.” 


For  example,  Figure  2.1  shows  the  B 2  beam  pattern  from  an  array  of  10 
receivers  spaced  0.5  wavelengths  apart.  These  beams  were  formed  using  Eq.  (2.25) 
with  the  received  pressure  waves  simulated  by 

Ri  = 


r  =  \l{X,ftf +  (ZSf 


(2.28) 


with  the  broadside  point  source  located  at  (0,  Zs )  for  Zs  equal  to  1,  10,  or  10,000 
wavelengths.  As  the  source  moves  closer  to  the  receiving  array,  the  assumption 
that  the  pressure  waves  are  planar  becomes  less  valid  and  eventually  the  beam- 
former  output  becomes  meaningless. 


Actually,  there  is  nothing  wrong  with  the  beamformer  —  we  merely  need 
to  modify  the  weights  of  Eq.  (2.26)  to  correspond  to  the  shape  of  the  pressure 
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wave  which  is  arriving  from  direction  9.  Furthermore,  in  some  cases  (such  as  the 
maximum  likelihood  processor),  it  would  be  desirable  to  have  the  beam  respond 
to  ail  pressure  waves  ai  riving  from  direction  0.  regardless  of  their  shape. 

An  adaptive  beamformer  could  be  constructed,  for  example,  which  re¬ 
sponds  to  a  spherically  spreading  pressure  wave  originating  in  the  direction  given 
by  9  at  some  distance  r.  The  key  component  of  such  a  processor  then  becomes  the 
method  used  to  estimate  the  range  r;  this  will  be  discussed  later  in  this  chapter. 
In  Chapter  4,  the  model  mismatch  problem  is  again  discussed  in  relation  to  the 
maximum  likelihood  processor. 

It  is  interesting  to  note  why  an  analogous  problem  does  not  arise  with 
the  holographic  approach.  As  with  a  beamformer,  Eq.  (2.8)  attempts  to  find  sine 
waves  across  an  equally  spaced  linear  array  (via  the  Fourier  transform.)  Any 
mismatch,  however,  is  perfectly  “undone”  by  the  inverse  Fourier  transform  which 
is  performed  after  the  propagation  effects  have  been  implemented  by  the  Green’s 
function. 

2.3  Receiver  Correlations 

Cross-correlations  between  the  signals  recorded  at  the  receivers  of  an  array 
can  be  used  to  locate  objects.  The  angle  at  which  the  object  lies,  relative  to  the 
array,  is  estimated  using  the  same  plane-wave  assumption  as  a  beamformer;  the 
angle  is  likewise  called  bearing.  However,  the  cross-correlation  method  can  in 
addition  estimate  range  by  using  the  curvature  of  the  wavefront.  Thus,  the  range 
should  be  great  enough  that  the  plane-wave  assumption  is  somewhat  correct  but 
not  so  great  that  there  is  no  detectable  wavefront  curvature. 
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Typically,  this  method  is  used  to  determine  the  location  of  objects  which 
are  the  source  of  the  pressure  wave;  however,  the  concept  can  be  adapted  to  find 
objects  which  only  scatter  the  wave  and  so  is  relevant  this  study. 

2.3.1  Basic  Correlation  Method 

This  derivation  will  use  the  geometry  shown  in  Figure  1.1  except  that 
there  shall  be  no  transmitter,  i.e.,  the  object  at  (Xs,  Zs)  is  itself  the  source  of 
the  pressure  wave.  To  simplify  the  equations,  we  shift  the  linear  array  so  that  the 
first  receiver  is  at  the  origin,  i.e.,  X R  =  0.  Cross  correlations  between  all  possible 
pairs  of  receivers  are  formed,  and  the  time  delay  diq  corresponding  to  the  highest 
peak  is  chosen  from  each  correlation.  The  subscripts  i  and  q  represent  any  two 
receivers.  This  relative  time  delay  can  be  written  in  terms  of  the  time  delay  at 
each  receiver  as  diq  =  d,  —  dq.  The  resulting  linear  system  of  equations  is  then 
solved  for  the  time  delays  d{  under  the  constraint  that  d\  =  0.  Thus,  we  can 
write  the  distance  r  between  the  first  receiver  and  the  object,  and  the  distance 
r,-  between  the  ith  receiver  and  the  object,  as 

r2  =  (X5)2  +  (ZS)\  (2.29) 

r?  =  (r  +  cdi)2  =  (Af  -  Xs)2  +  (Zs)2  ,  (2.30) 

respectively.  We  next  subtract  two  instances  of  the  latter  equation,  for  receivers 
i  and  q,  after  first  multiplying  each  by  the  receiver  coordinate  of  the  other.  This 
produces 


R  -2  _  XR  2 

*■  1  '  n 


Xqrt 


(2.31) 


=  XqR  r 2  +  IrXfcdi  +  XRc2d] 


XR  r2  -  2rA 


R 


cdn  — 


X,Rc2d2 
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=  Xf  (xf )2  -  2XsXRXf  +  X*  (X5)2 
-  XR  (XR)2  +  2XsXRXf  -  Xf  (Xsy 

+  (xqR-xf)  (zs)2 


2 r  (XRcdi  -  Xf cdq)  +  (xfc2d2  -  XRc2d2) 

=  xfxf  (xf  -  Xf)  +  (Xf  -  Xf)  (r2  -  (A-)2  -  (Zs)2)  . 


(2.32) 


Eq.  (2.29)  can  be  used  to  eliminate  the  last  term  in  Eq.  (2.32)  so  that  we  can 


solve  for  r: 


X?X?  (Af  -  X«)  -  c2  (A'ffi  -  A :£4) 
2c(X*d,-X!>d,) 


(2.33) 


When  there  are  only  three  equally  spaced  receivers  with  spacing  L,  we  can  let 
Xf  =  L  and  Xf  =  —  L  to  yield 


L2  -  (<?  +  rf2) 

c(di  +  dq) 


(2.34) 


In  order  to  use  all  available  information,  we  average  the  estimate  of  r  from 
Eq.  (2.33)  over  all  possible  pairs  to  obtain 


wf  Af  A*  (A/2  -  A'*)  -  c2  (A*4  -  A 

-  h  h  ^  (*,**  -  *,%) 


(2.35) 


where  normalization  is  accomplished  by  letting  Wiq  =  1/(nr  ( NR  —  l)Y 


If  we  can  assume  that  the  received  waves  are  planar,  the  bearing  angle  0 
as  shown  in  Figure  1.1  can  be  computed  from  each  time  delay  and  the  results 


(2.36) 


The  optimum  estimation  of  range  and  bearing  using  correlations  was  ad¬ 
dressed  by  Hahn27  for  linear  arrays  of  receivers.  The  object  producing  the  pres¬ 
sure  wave  was  assumed  to  emit  wide-band,  zero-mean  Gaussian  noise,  while  ad¬ 
ditional,  independent  noise  was  introduced  at  each  receiver. 

In  Hahn’s  results,  the  weights  are  chosen  to  not  only  provide  normalization 
but  to  also  minimize  the  error  covariance  matrix  and  thus  approach  the  “opti¬ 
mum”  Cramer-Rao  bound,  the  theoretically  optimum  performance  of  an  array. 
Hahn  states  that  the  calculation  of  the  weights  is  trivial  in  the  numerical  case 
but  complicated  for  most  theoretical  cases.  He  provides  a  formula  for  the  weights 
only  for  the  case  of  arrays  arranged  symmetrically  around  the  origin  where  the 
noise  level  is  the  same  at  each  receiver. 

2.3.3  Comments  on  the  Correlation  Method 

The  ambiguity  function  corresponding  to  each  time  delay  di  used  in  the 
correlation  method  is  a  straight  line  given  by  one  term  of  Eq.  (2.36).  The  range 
estimate  of  Eq.  (2.35)  can  therefore  be  viewed  as  a  calculation  of  where  two  of 
these  lines,  as  defined  by  di  and  d?,  intersect. 


2S 

averaged  to  yield 


Nr  j.  \ 

0  =  arcsin  [  Y, 

i=l  Ai  J 


where  normalization  is  accomplished  by  letting  —  1  /NR. 


2.3.2  Optimum  Correlation  Method 


From  this  viewpoint  it  is  easy  to  see  where  failure  can  occur:  if  the  two  time 
delays  used  in  each  term  of  Eq.  (2.35)  are  incompatible,  the  resulting  two  lines 
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may  cross  at  a  location  far  removed  from  the  actual  object.  The  incompatibility 
can  occur  because  the  delays  diq  are  calculated  by  picking  the  highest  peak  in  the 
correlation;  this  picking  process  may  not  be  robust  when  there  are  several  nearly 
equal  peaks  such  as  would  occur  when  there  are  several  objects  or  an  extended 
object  emitting  pressure  waves.  That  is,  even  if  dt  is  due  to  one  object  while  dq  is 
due  to  another,  the  algorithm  will  continue  attempt  to  determine  a  single  range 
and  a  single  bearing  in  each  term  of  Eq.  (2.35). 

Hahn  notes  that,  if  more  than  one  source  exists,  the  signals  must  be 
preprocessed  so  as  to  introduce  nulls  in  the  direction  of  the  additional  sources. 
This  is  perhaps  the  greatest  weakness  of  this  approach  and  is  a  strong  motivation 
for  seeking  other  methods  which  use  similar  correlation  techniques  but  which  can 
better  accommodate  multiple  or  extended  objects. 

2.4  Pattern  Matching 

Even  before  one  begins  to  consider  the  concepts  of  super  resolution,  it 
becomes  apparent  that  the  traditional  methods  of  the  previous  sections  have 
difficulty  accommodating  additional  features  such  as  a  priori  information  about 
the  objects,  multiple  frequencies,  or  more  complex  propagation  models.  It  may 
be  that  these  added  features  could  not  have  been  seriously  considered  before  the 
arrival  of  the  modern  digital  computer,  since  they  would  have  been  difficult  or 
impossible  to  implement  with  simple  sum-and-delay  beamformers,  or  with  the 
optical  reconstruction  methods  used  in  the  early  days  of  holographic  imaging. 

One  has  the  sense  that  the  added  features,  as  well  as  the  traditional  meth¬ 
ods,  should  all  fit  into  some  general  model  since  they  are  all  based  upon  the  wave 
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equation’s  ability  to  describe  the  propagation  of  pressure  waves  through  a  ho¬ 
mogeneous  acoustic  medium.  In  the  following  sections,  such  a  general  model  is 
constructed  around  the  concept  of  a  fictitious  scatterer;  then,  several  special  cases 
of  the  general  model  are  considered,  two  of  which  correspond  to  the  traditional 
holographic  and  beamforming  methods. 

2.4.1  Fictitious  Scatterers 

It  is  possible  that  several  workers  have  proposed  the  approach  discussed 
in  this  section;  however,  the  author’s  introduction  to  this  approach  came  from  a 
particularly  clear  paper  published  by  Ermert  and  Karg28  which  described  their 
theoretical  and  experimental  development  of  what  they  called  multifrequency 
acoustical  holography. 

They  begin  by  noting  that  the  theory  for  a  matched  filter  is  the  basis  for 
their  approach.  First,  they  assume  that  only  a  single  scatterer  at  some  specific 
location  is  present  and  then  calculate  what  this  fictitious  scatterer  would  have 
yielded  at  their  receivers.  Second,  they  correlate  the  complex  conjugate  of  the 
signal  from  this  fictitious  scatterer  with  the  signals  actually  received.  And  third, 
the  fictitious  scatterer  is  moved  through  all  locations  in  object  space,  i.e.,  all  lo¬ 
cations  which  could  possibly  contain  an  actual  scatterer.  The  locations  at  which 
we  find  peaks  in  the  output  of  the  correlation,  then,  are  interpreted  as  the  loca¬ 
tions  of  the  true  object  or  objects.  The  correlation  \k,  from  their  Eq.  (8),  would 
be  expressed  in  the  notation  of  this  study  for  continuous  time  and  a  continuous 
receiving  aperture  as 

+oo  -1-00 

J  J  RH  (x,t)  •  (i?F(x,f)) 

— oo  — oo 


^(XS,ZS)  = 


dx  dt 


(2.37) 


where  RH  and  RF  are  the  actual  and  fictitious  signals,  respectively,  at  receiver 
aperture  locations  XR .  The  vertical  bars  (| . . .  |)  denote  the  absolute  value  op¬ 
eration.  The  fictitious  scattering  object  is  at  (Xs,  Zs).  When  the  aperture  and 
time  are  sampled,  the  infinite  integrals  are  replaced  by  finite  sums  to  yield 

Nr  Nt 

1(A'S,ZS)=  AX  At  ZZRlt  (Rf,)’  .  (2.38) 

1=1  /=1 

Through  the  concept  of  a  fictitious  scatterer,  the  power  and  generality  of 
matched  filter  theory,  originally  developed  for  analyzing  time  histories,  has  been 
extended  to  analyze  the  space-time  patterns  received  at  the  array.  Virtually  any 
type  of  transmitted  signal  can  be  accommodated  (Ermert  and  Karg  used  a  linear 
frequency  sweep)  as  well  as  any  spatial  arrangement  of  receivers. 

Furthermore,  any  sequence  of  events  that  affect  the  propagation  of  the 
pressure  wave  can  be  incorporated.  The  sequence  of  events  can  include  closed- 
form  formulas;  for  example,  RF  could  be  calculated  for  simple  free-space  prop¬ 
agation  using  Eq.  (1.2).  But,  more  importantly,  sequences  of  events  which  can, 
say,  be  described  only  numerically  can  also  be  incorporated.  Even  propagation 
effects  known  only  empirically  due  to,  say,  preliminary  system  tests29,30  could 
be  used  to  generate  the  fictitious  receiver  signals  for  the  desired  locations  of  the 
fictitious  scatterers. 

Some  of  these  variations  will  be  examined  in  subsequent  sections  of  this 


chapter  and  in  Chapter  5. 
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2.4.2  Reconstruction  as  Pattern  Matching 

When  a  single  frequency  sine  wave  is  used,  the  integral  over  time  in 
Eq.  (2.37)  merely  multiplies  the  result  by  a  constant  that  can  be  dropped.  The 
pattern-matching  formula  thus  becomes 


T(A5,Z5)  = 


+  0O 


J  Rh{x)-  (RF(x)Y  dx 


(2.39) 


The  reconstruction  represented  by  Eq.  (2.8)  can  be  rewritten  as  an  inverse  Fourier 
transform  or  a  convolution: 

+oo 

R(XS ,  Zs)  =  — y==  J  Rh (kx)  ■  exp  (j  kz  Zs  —  j  kx  i)  dkx,  (2.40) 

— oo 

=  (. RH(kx )  ■  exp  (j  kz  ZS)) ,  (2.41) 

+0O 

=  J  RH(x)-G(Xs -x,Zs)  dx  (2.42) 


where  fyG  (x,  Z5)  j  =  exp  (j  kz  Zs In  this  comparison  between  holography 
and  pattern  matching,  it  is  sufficient  to  consider  only  the  non-evanescent  waves; 
for  this  case,  Cohen31  has  shown  that  G  for  a  given  Zs  (or,  in  his  notation,  G* 
for  a  given  D),  is  equal  to 


(2.43) 


where  r2  =  x2  - f  {Zs^  .  When  r  »  A,  terms  on  the  order  of  p  are  negligible 
compared  to  terms  on  the  order  of  j.  Thus, 


(2.44) 


Upon  examining  Eq.  (2.39)  and  Eq.  (2.42)  we  see  that 


G(xs  ~x,Zs)  =  (RF(x)Y  (2.45) 

or 

Rr(x)  =  ->jL(Zl)  eM-jkr) 

2  K  \  r  J  r 

2  2 

where  x  — ►  Xs  —  x  so  that  we  now  have  r~  =  (X5  —  x)  4-  (-Z5)  .  A  comparison 
of  Eq.  (2.39)  and  Eq.  (2.42)  indicates  that  the  two  formula  are  identical  when  the 
fictitious  source  is  given  by  Eq.  (2.46).  Although  not  explicitly  used  in  Eq.  (2.8), 
an  absolute  value  operation  is  typically  added  when  the  image  is  plotted. 

The  source  radiation  expressed  by  Eq.  (2.46)  consists  of  spherically  spread¬ 
ing  waves  with  a  radiation  pattern  given  by 

Zs 

cos  (</>)  =  —  (2-47) 

where  (f>  =  0  points  towards  the  negative  Z  axis.  As  explained  by  Skudrzyk,32 
this  pattern  is  due  to  the  fact  that  for  the  development  of  the  Helmholtz-Huygens 
integral,  the  equivalent  sources  at  the  hologram  needed  to  have  zero  radiation  in 
the  backwards  direction.  This  was  accomplished  by  equating  each  of  these  sources 
to  the  sum  of  a  point  source  and  a  dipole,  with  the  consequence  that  a  cardioid 
radiation  pattern  was  generated. 

In  some  applications  of  holography,  such  as  nearfield  studies,  the  angle  <f> 
between  a  source  of  acoustic  energy  and  a  point  in  the  hologram  may  become 
quite  large,  or  even  approach  90  degrees.  In  these  cases,  the  fictitious  source 
used  by  the  holographic  equation,  with  its  cardioid  radiation  pattern,  may  do  a 
poor  job  of  modeling  the  more  omnidirectional  radiation  from,  say,  a  point  on  a 
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vibrating  surface.  Thus,  the  pattern-matching  approach  could  provide  a  distinct 
advantage  by  allowing  us  to  use  an  omnidirectional  source  in  the  analysis. 


The  two  methods  share  the  ambiguities  in  the  interpretation  of  their  im¬ 
age.  Either  can  be  used  to  estimate  the  object  field  at  arbitrary  points 
although  holography  can  product  a  more  efficient  calculation  when  the  object 

field  is  evaluated  at  the  same  locations  along  the  X  axis  as  the  receivers,  i.e.,  at 
yS  __  yR 

Am  —  Aj-  . 


As  shown  by  Ermert  and  Karg,  their  integral  over  time,  combined  with  a 
sufficiently  large  bandwidth,  can  reduce  the  range  ambiguity  inherent  in  holog¬ 
raphy.  This  capability  will  be  examined  again  in  Chapter  5. 


2.4.3  Beamforming  as  Pattern  Matching 

A  comparison  of  the  Bartlett  beamformer,  as  given  by  Eq.  (2.25),  and  the 
pattern  matching,  as  given  by  Eq.  (2.39),  shows  that  the  fictitious  source  used  in 
beamforming  is  given  by 


R[  =  exp  j  kXf1  sin(0)j 


(2.4S) 


which  is  of  course  just  the  pressure  which  would  have  been  received  from  a  source 
radiating  plane  waves. 


2.4.4  Correlation  Method  as  Pattern  Matching 

The  pattern-matching  approach  can  be  categorized  as  a  hypothesis  tester: 
one  selects  a  suitable  domain  of  possible  object  locations,  and  then  tests  each 
hypothesis  with  Eq.  (2.38).  On  the  other  hand,  the  correlation  method  provides 
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an  immediate  estimate  of  the  object  location.  The  advantage  in  such  a  direct 
estimate  is  that  fewer  calculations  are  needed  to  reach  a  conclusion;  the  disad¬ 
vantage  is  that  it  can  fail  disastrously  when  there  is  more  than  one  object,  i.e., 
there  exists  no  single  answer.  In  these  situations,  the  pattern-matching  approach 
would  seem  to  be  more  robust  as  it  can  handle  multiple  sources. 

Because  of  the  fundamental  difference  in  the  approach  used  by  the  two 
methods,  as  noted  in  the  previous  paragraph,  forming  a  direct  mathematical  com¬ 
parison  may  be  impossible.  An  initial  comparison  of  the  two  methods  would  sug¬ 
gest  that  the  integral  of  the  received  and  fictitious  patterns  over  time  is  the  main 
source  of  range  information,  and  as  such  is  somewhat  analogous  to  the  correlation 
method’s  calculation  of  range,  Eq.  (2.35).  The  integral  over  space  is  somewhat 
analogous  to  the  correlation  method’s  calculation  of  bearing,  Eq.  (2.36). 

A  further  comparison  of  how  the  time-history  information  is  used  can 
be  formed  by  first  simplifying  the  algorithm:  assume  that  the  bearing  0  =  0; 
and  that  changes  in  the  amplitude  of  the  received  signals  across  the  array  are 
negligible.  This  means  that  we  only  need  to  seek  the  value  of  the  unknown  range 
r,  which  is  equivalent  to  seeking  the  value  of  the  time  delays  since  the  speed  of 
sound  is  a  known  constant. 

When  these  assumptions  are  applied  to  the  correlation  method,  Eq.  (2.35), 
we  find  that  the  information  obtained  through  the  inter-element  time  delays 
diq  (the  location  of  the  peak  in  the  cross-correlation  between  the  signals  from 
receivers  i  and  q )  can  now  be  obtained  instead  through  the  time  delays  dn  (the 
location  of  the  peak  in  the  cross-correlation  between  the  signals  from  receivers 
i  and  1).  This  also  eliminates  the  need  to  solve  the  simultaneous  system  of 
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equations  since  dt\  =  di~d\  and  d\  =  0.  We  can  roughly  summarize  the  resulting 
process  as  follows  where  the  first  two  steps  constitute  the  cross  correlation: 

(1)  For  some  time  delay  r,-,  evaluate  the  match 
using  £/= i  Ri  ( U )  ■  R*  ( ti  +  r,); 

(2)  Repeat  for  all  time  delays  where  Tm;n  <  Ti  <  Tmax. 

(3)  Let  d{  =  T{  —  t\  using  the  value  of  Ti  that  produced  the  best  match; 

(4)  Repeat  for  all  receivers  i  where  1  <  i  <  NR;  then 

(5)  Calculate  the  range  by  inserting  the  d,-  into  Eq.  (2.35). 

Turning  now  to  the  pattern-matching  method,  we  can  roughly  summarize 
its  algorithm  as  follows: 

(1)  For  some  range  rtest ,  calculate  the  corresponding  fictitious  r,-  at 
receiver  i; 

(2)  Repeat  for  all  receivers  i  where  1  <  i  <  NR-, 

(3)  Evaluate  the  match  using  Eq.  (2.38); 

(4)  Repeat  for  all  ranges  where  rmin  <  rtest  <  ^moi!  then 

(5)  Let  r  equal  the  value  of  rtest  that  produced  the  best  match. 

Both  methods  use  two  loops:  one  over  receivers,  and  one  over  either  time 
delays  (step  2  in  the  correlation  method)  or,  equivalently,  over  range  (step  4  in 
pattern  matching).  The  result  is  then  merely  the  value  that  produced  the  best 
match.  The  main  difference  between  the  two  methods,  therefore,  seems  to  be 
that  the  pattern-matching  method  has  reversed  the  order  of  its  two  loops  relative 
to  the  correlation  method. 
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2.5  Array  Geometry 

In  this  study,  we  have  limited  our  investigations  to  processing  methods 
which  use  a  linear  array.  However,  it  seems  prudent  to  at  least  review  other 
methods,  using  two-  or  three-dimensional  arrays,  to  determine  if  there  might 
be  some  techniques  or  beneficial  receiver  spacing  which  could  be  applied  to  a 
one-dimensional  array. 

An  explicit  inversion  of  the  Helmholtz  equation  was  used  by  Ball  et  al ,33 
to  develop  a  method  for  reconstruction  using  measurements  chosen  to  lie  on  the 
surface  of  a  sphere.  The  resulting  expressions  were  in  terms  of  spherical  harmonics 
and  could  accommodate  both  plane  wave  and  spherical  wave  ultrasonic  sources. 
Norton  and  Linzer34  investigated  the  exact  inverse  scattering  solution  for  plane, 
cylindrical,  and  spherical  arrays.  However,  in  order  to  simplify  the  mathematics, 
they  only  considered  backscattering;  i.e.,  each  receiver  individually  operated  as 
first  a  transmitter  and  then  a  receiver  so  that  a  given  receiver  never  encountered 
waves  arising  from  the  other  receiver-transmitter  locations.  In  both  of  these 
studies,  the  goal  was  to  estimate  the  velocity  and  hence  the  density  of  the  media 
within  the  spherical  array.  In  neither  case  were  any  particular  requirements  for 
the  spacing  or  arrangement  of  the  receivers  in  the  array  noted. 

The  reconstruction  of  a  three-dimensional  field  from  a  partial  sphere  of 
measurements  with  a  low  signal-to-noise  ratio  was  investigated  by  Bresler  and 
Macovski.35  Their  statistical  description  of  the  object  of  interest  and  its  rela¬ 
tionship  to  the  measurement  locations  has  the  potential  of  offering  some  sound 
guidelines  to  array  configurations;  unfortunately,  they  did  not  address  that  issue. 
The  reconstruction  of  three-dimensional  domains  using  two-dimensional  arrays 


3S 


was  also  considered  by  Koppelmann  and  Keating36  for  a  rectangular  array  and 
by  Norton37  for  a  circular  array. 

A  review  of  two-dimensional  arrays  by  Nigam38  provides  insights  which 
are,  after  15  years,  still  relevant.  The  number  and  size  of  the  elements,  the 
overall  size  of  the  array,  and  the  electronic  challenge  of  connecting  large  numbers 
of  transducers  were  discussed;  however,  only  fully  populated  arrays  with  equally 
spaced  elements  were  considered. 

One  of  the  most  comprehensive  studies  of  a  linear  array  was  presented  by 
Carter.39  In  particular,  he  examined  the  case  wherein  the  array  was  constrained 
to  be  of  a  given  length  L  with  a  given  number  of  receivers  M.  Ideally,  M  is 
large  enough  so  that  there  is  a  receiver  every  half- wavelength.  Carter  considered 
the  question  of  how  to  optimally  distribute  the  receivers  when  less  than  the  ideal 
number  were  available.  When  the  goal  was  to  minimize  the  variance  of  the  bearing 
estimate,  the  optimum  configuration  was  shown  to  be  one  where  half  the  receivers 
were  at  each  end  of  the  array.  If  instead  the  goal  is  to  minimize  the  variance  of 
the  range  estimate,  the  optimum  arrangement  was  shown  to  be  one  with  half 
the  elements  at  the  center  with  one-quarter  of  the  elements  at  each  end.  To 
minimize  the  area  of  uncertainty  (roughly,  the  product  of  the  range  and  bearing 
variance),  Carter  found  numerically  that  the  optimum  arrangement  was  one  with 
one-third  of  the  receivers  in  the  center  and  one-third  at  each  end.  However,  at 
very  short  ranges,  the  minimum  area  of  uncertainty  was  shown  to  result  using 
the  arrangement  that  minimized  the  bearing  variance  while  for  very  large  ranges 
the  arrangement  that  minimized  the  range  variance  was  recommended.  In  every 
case,  the  receivers  were  spaced  a  half-wavelength  apart  in  each  subarray. 
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A  similar  question  was  addressed  by  Pillai  et  al.40  who  sought  the  optimum 
placement  of  receivers  in  a  linear  array  for  use  in  spatial  spectrum  estimation  (i.e., 
bearing  estimation)  using  one  of  the  methods  based  upon  the  eigenstructure  of 
the  covariance  matrix  (discussed  further  in  the  next  chapter).  It  was  shown 
that  the  M  receivers  should  be  placed  at  locations  equal  to  a  half- wavelength 
multiplied  by  a  set  of  integers  known  as  a  Caratheodory  sequence.  For  example, 
when  M  =  4,  the  receivers  should  be  located  at  0,  £A,  |A,  and  £A. 

Two  conclusions  seem  to  be  in  order.  First,  the  optimum  array  geometry 
and  receiver  spacing  is  seldom  a  universal  one;  it  instead  depends  upon  which 
parameter  one  wishes  to  optimize  and  which  processing  method  is  to  be  used 
with  the  resulting  array.  Second,  when  using  a  one-dimensional  linear  array, 
there  seems  to  be  no  method  which  claims  to  improve  the  results  beyond  what 
can  be  obtained  when  the  array  is  fully  populated  from  end  to  end  with  elements 
spaced  every  half-wavelength.  Since  the  simulations  used  in  Chapters  4  and  5 
always  use  such  an  ideal  array,  the  issue  of  receiver  spacing  does  not  arise. 

2.6  Time-of-Flight  Information 

One  strong  motivation  for  not  using  the  full  time-history  of  the  signal  is 
that  it  is  difficult  to  analytically  represent  processing  operations  on  such  a  signal. 
Nevertheless,  it  illustrative  to  consider  what  is  being  lost  by  not  using  the  full 
time-history. 

The  simplest  system  for  determining  the  location  of  a  scattering  object  is 
a  single  transmitter  and  a  single  receiver.  The  general  expression  for  given  in 


Eq.  (1.2),  becomes 


Ri(ti)  =  Sx 


E(t-  Tf) 
D\  •  Dt 


(2.49) 


where  S\  is  the  reflection  coefficient.  The  time- of- flight  (TOF)  Tf  and  the  dis¬ 
tances  were  defined  in  Section  1.4.  The  amplitude  of  R\  is  relatively  insensitive 
to  small  changes  in  the  distances,  and  includes  the  unknown  coefficient  S\.  The 
TOF  is  related  to  the  length  of  the  path  traveled  by  the  pressure  wave  through 


F  D\  +  Dt 


Tf  = 


(2.50) 


Unfortunately,  this  value  is  identical  for  all  points  on  an  ellipse  which  has  its  foci 
at  the  transmitter  and  receiver.  That  is,  in  an  experiment,  we  do  not  know  D\ 
and  Dt  separately  —  we  only  know  their  sum  by  using  the  measured  TOF  and 
Eq.  (2.50).  Thus,  knowing  the  TOF  is  not  the  same  as  knowing  the  range  of  the 
scatterer  since  the  range  or  distance  is  different  to  each  point  on  the  ellipse. 


The  next  logical  step  is  to  add  either  additional  transmitters  or  addi¬ 
tional  receivers,  where  one  ellipse  of  ambiguity  is  contributed  by  each  transmitter- 
receiver  pair.  If  additional  transmitters  are  used,  we  must  be  able,  at  each  re¬ 
ceiver,  to  determine  which  transmitter  was  the  source  in  order  to  calculate  the 
ellipse  of  ambiguity  from  the  time  delay.  This  can  be  done  either  by  encoding 
the  signals,  or  by  activating  the  transmitters  one  at  at  time. 


A  less  complicated  system  can  be  formed  by  instead  using  only  one  trans¬ 
mitter  and  multiple  receivers.  In  general,  all  ellipses  will  add  coherently  only 
at  the  one  location  which  corresponds  to  the  scatterer.  However,  there  will  in 
general  be  many  secondary  peaks  generated  where  any  two  ellipses  cross.  When 
there  are  multiple  scattering  objects,  each  contributes  its  own  set  of  ellipses  of 
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ambiguity  to  the  image. 

The  previous  paragraphs  have  assumed  that  th<  POF  has  been  used. 
Many  imaging  methods,  however,  do  not  use  this  time  directly;  instead,  they 
assume  that  the  pressure  waves  consist  of  a  single  frequency  and  then  represent, 
in  effect,  the  TOF  as  the  phase  of  the  complex  signal 


Ri  =  Si 


k  (D\  +  Dt )) 
D\  •  Dt 


(2.51) 


While  this  is  fundamentally  the  same  information,  a  phase  ambiguity  of  2n  has 
been  introduced.  This  is  equivalent  to  a  time  ambiguity  of  one  period  or  a 
distance  ambiguity  of  one  wavelength  (A).  Thus,  instead  of  one  ellipse  for  each 
transmitter-receiver  pair,  there  are  an  infinity  of  concentric  ellipses  where  the 
TOF  for  any  two  adjacent  ellipses  differ  by  one  period. 


Some  methods  are  built  upon  a  correlation  between  the  signal  at  each 
receiver  and  the  transmitted  signal.  A  full  cross-correlation,  using  a  full  set  of 
inner-signal  delays,  would  identify  the  TOF  as  the  value  of  the  inner-signal  delay 
at  which  the  best  correlation  was  obtained.  However,  in  most  cases,  the  method 
uses  only  the  zero-delay  correlation  —  that  is,  they  only  examine  how  well  the 
two  signals  being  compared  match  “as  is”  by  calculating  only  one  correlation 
with  an  inner-signal  delay  of  zero.  Thus,  ambiguities  are  again  introduced.  A 
similar  ambiguity  also  results  when  the  correlations  are  between  receivers.  And, 
as  noted  in  Section  2.3.3,  the  process  of  choosing  the  one  peak  which  yields  the 
TOF  may  not  be  robust. 


The  conclusion,  then,  is  that  any  processing  method  that  reduces  the 
TOF  information  from  the  receivers  to  a  single  complex  number  is  discarding 
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information  that  could,  potentially,  be  used  to  improve  the  resolution.  Although 
the  simulations  used  in  Chapters  4  and  5  do  not  use  this  TOF  information,  a 
brief  discussion  of  how  to  incorporate  the  TOF  into  the  matching  methods  will 
be  given  in  Chapter  6. 
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Chapter  3 

High  Resolution  Methods 

The  methods  for  processing  signals  from  an  array  which  are  examined  in 
this  chapter  will  include  ones  seeking  high  resolution,  or  the  ability  to  distinguish 
two  closely  space  objects,  as  well  as  methods  which  seek  high  accuracy,  or  the 
ability  to  determine  the  location  of  a  single  object  with  high  precision.  This 
type  of  accuracy  is  sometimes  also  described  as  the  variability  of  the  location  or 
the  variance  of  the  range  and  bearing  estimates.  Although,  technically  speaking, 
these  are  different  goals,  the  methods  for  achieving  them  are  similar  enough  that 
we  can  refer  to  all  the  methods  as  high  resolution  ones. 

3.1  Spectral  Estimators 

In  the  analysis  of  the  spectrum  of  a  time  series  —  that  is,  a  sequence  of 
samples  of  some  phenomena  versus  time  —  the  simplest  method  for  obtaining  the 
spectrum  is  to  form  its  Fourier  transform.  This  results  in  a  spectral  resolution 
of  l/T  where  T  is  the  length  of  the  sampling  period.  A  high-resolution  spectral 
estimator  attempts  to  generate  a  result  which  has  the  same  effect  as  an  increase 
in  T  with  a  corresponding  increase  in  the  resolution  of  the  spectrum. 

This  idea  can  be  transferred  to  line  arrays  by  letting  the  location  along 
the  array  take  on  the  role  of  time  and  the  vector  wavenumber  take  on  the  role  of 
frequency.  Then,  the  enhancements  generate  a  result  which  has  the  same  effect  as 
an  increase  in  the  array  length  L  with  a  corresponding  increase  in  the  resolution 
of  the  spatial  spectrum.  Thus,  any  spectral  estimator  that  can  be  used  to  enhance 
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the  temporal  spectrum  can  be  applied  to  the  spatial  spectrum.  Several  of  the 
methods  discussed  later  in  this  chapter  can  be  classified  as  spectral  estimators, 
and  were  originally  developed  for  the  time  domain. 

It  is  important  to  not  confuse  this  method,  which  increases  the  resolution 
of  the  spectrum,  with  those  that  instead  extend  the  spectrum  to  higher  spatial 
frequencies.  The  latter  effect  is  the  same  as  taking  the  samples  closer  together  in 
space.  Such  an  effect  cannot  yield  any  additional  information  when  the  samples 
are  already  half  a  wavelength  or  closer  together  except  when  evanescent  waves 
are  involved. 

It  is  also  important  to  note  that  most  spectral  estimators,  including  the 
FFT,  produce  estimates  which  are  equally  spaced  in  the  wavenumber  domain  but 
not  in  angle.  For  example,  in  the  two-dimensional  space  used  in  this  study,  we 
would  have  estimates  at  kx  =  0,  ±Afcz,  ±2 A kx,...  where  A kx  —  2ir / L  for  an 
array  of  length  L.  Using  the  angular  spectrum  interpretation  of  Goodman,14  each 
of  these  wavenumbers  correspond  to  a  plane  wave  of  frequency  /  and  wavenumber 
k  arriving  at  angle  0  so  that  kx  =  k  sin(0)  or  0  =  arcsin (kx/k).  Thus,  the 
estimates  are  relatively  far  apart  around  0  —  0  and  most  closely  spaced  around 
0  =  ±90°.  Unfortunately,  we  often  are  seeking  the  highest  accuracy  in  bearing 
around  broadside  or  0  =  0. 

3.1.1  ARMA  Methods 

One  of  the  more  widely  used  class  of  spectral  estimators  consist  of  either  a 
polynomial,  the  inverse  of  a  polynomial,  or  the  ratio  of  two  polynomials.  Equiva¬ 
lent  names  for  the  polynomial  estimator  are  all-zero,  moving  average  (MA),  and 
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feedforward  filter  while  equivalent  names  for  the  inverse-polynomial  estimator  are 
all-pole,  autoregressive  (AR),  and  feedback  filter.  An  ARMA  estimator  combines 
both  and  is  the  ratio  of  two  polynomials.  The  all-zero  and  all-pole  designations 
refer  to  the  Laplace  transform  of  the  filter.  The  coefficients  of  these  polynomials 
are  first  determined  using  the  available  data;  then,  the  polynomials  can  be  used 
to  generate  additional  synthetic  data  between  existing  data  points  (interpolation) 
or  beyond  the  domain  of  the  existing  data  points  (extrapolation). 

In  addition  to  the  choice  to  use  an  AR,  MA,  or  ARMA  model,  one  must 
also  select  the  method  for  determining  the  order  of  the  polynomial(s),  the  method 
for  the  calculation  of  the  polynomial  coefficients  (either  directly  or  through  iter¬ 
ation),  and  the  subject  of  estimation. 

3.1.2  Options  for  the  Subject  of  Estimation 

One  of  the  difficulties  in  assembling  this  review  of  ARMA  methods  has 
been  the  fact  that  the  polynomials  can  be  fitted  to,  and  used  to  estimate,  several 
different  aspects  of  the  signal.  The  most  useful  basis  on  which  to  categorize  these 
variations  seems  to  be  whether  the  estimator  operates  on  time-domain  data  (i.e., 
the  time-history  of  the  signal  at  each  receiver)  or  on  space-domain  data.  In  the 
former  case,  the  improvements  are  equivalent  to  sampling  the  signal  either  more 
rapidly  or  for  a  longer  period  of  time.  Since  neither  of  these  directly  affect  the 
processing  of  data  from  an  array,  they  shall  not  be  considered  further. 

Estimators  which  operate  on  space-domain  data  require  that  the  informa¬ 
tion  at  each  receiver  first  be  represented  as  a  complex  number  —  that  is,  the 
time-history  of  the  signal  cannot  be  used  directly.  Thus,  the  ambiguities  dis- 
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cussed  in  Section  2.6  are  inherent  in  ARMA  methods.  These  estimators  can  also 
operate  either  directly  on  the  data,  or  on  some  function  of  the  data.  For  exam¬ 
ple,  some  fit  the  polynomial  'H(i)  to  the  ratio  of  the  complex  signal  at  adjacent 
receivers  R,  and  Ri+i  so  that 


Rj+x 

Ri 


=  y(k) 


for 


i  =  1,2,. ..,1V. 


(3.1) 


Other  versions  may  operate  on  the  Fourier  transform  of  the  signals,  i.e.,  in  the 
wavenumber  domain. 


The  literature  which  has  been  reviewed  suggests  that  there  are  several 
goals  to  be  considered  when  choosing  the  data  or  function  of  the  data  which  is 
to  be  the  subject  of  estimation  (SE): 

•  Some  technique  must  exist  for  the  selection  of  the  correct  model  (MA, 
AR,  or  ARMA)  for  the  SE. 

•  Some  technique  must  exist  for  the  selection  of  the  correct  order  for  the 
polynomial(s)  to  be  fitted  to  the  SE. 

•  There  must  be  a  closed-form  or  convergent  iterative  method  for  evaluating 
the  coefficients  from  existing  SE  data. 

•  The  SE  should  be  resistant  to  noise  and  measurement  errors. 

As  an  example  of  this  last  goal,41  it  has  been  found  that  the  average  of  the  cross¬ 
correlation  between  receiver  pairs  is  more  robust  than  the  cross-correlation  of  the 
average  receiver  signal  since,  in  the  former,  the  cross-correlation  helps  remove 
noise  before  the  averaging  process  is  applied.  The  optimum  compromise  between 
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these  goals  is  not  always  obvious. 

3.1.3  Choosing  the  Estimator  Parameters 

The  error  in  the  estimation  will  be  minimized  when  the  model  and  order 
match  the  underlying  physical  process.  A  study  of  11  categories  of  methods  for 
the  estimation  of  the  power  spectral  density  of  discrete  time  series  by  Kay  and 
Mar  pie4'  concluded  that  "if  the  model  is  inappropriate  ...  poor  (biased)  spectral 
estimates  will  result.”  The  question  of  interest  is  whether,  in  the  absence  of  any 
prior  knowledge,  the  model  and  order  can  be  derived  robustly  from  the  data. 

The  MA  and  ARMA  methods  were  used  by  Abdel-Aal  et  al ,43  to  estimate, 
in  the  space  domain,  the  ratio  of  the  signal  from  adjacent  receivers.  The  source 
of  the  acoustic  waves  in  their  computer  simulations  was  either  one  or  two  point 
sources.  They  found  that  the  MA  model  led  to  errors  so  large  that  it  was  not 
useful  when  there  was  more  than  one  source.  While  the  ARMA  model  yielded 
better  results,  no  discussion  on  the  choice  of  model  order  is  provided.  They 
seemed  to  have  simply  tried  different  model  orders  until  the  predictions  matched 
what  they  knew  beforehand. 

The  MA,  AR  and  ARMA  models  were  used  by  Gutowski  et  a/.44  to  predict 
the  spectrum  of  three  time  series  which  were  designed  to  be  purely  MA,  purely 
AR,  or  ARMA.  Their  simulations  revealed  that  using  a  model  which  is  incorrect 
for  the  data  can  result  in  predicted  spectra  which  have  the  wrong  shape  or  peaks 
which  do  not  coincide  with  the  true  spectral  peaks.  In  conclusion,  they  state: 

One  problem  which,  in  our  opinion,  remains  essentially  unsolved,  is 

a  practical  means  to  determine  a  priori  whether  a  real  life  situation 


corresponds  to  an  AR,  MA,  or  ARMA  process.  Another  problem 
is  how  to  effectively  determine  the  order  of  a  given  process  other 
than  (by)  the  empirical  methods  used  here.  For  the  AR  model, 
the  order  of  the  feedback  component  can  be  determined  from  the 
behavior  of  the  partial  correlation  coefficient.  This  test,  however, 
often  breaks  down  when  applied  to  the  ARMA  model. 

More  recently,  Wax45  noted  that  two  of  the  methods  which  have  been 
proposed  for  estimating  the  order  of  an  AR  model  (called  AIC  and  MDL)  are 
appropriate  only  when  there  are  a  large  number  of  data  samples  (because  the 
estimates  are  asymptotic)  and  the  data  can  be  modeled  as  Gaussian.  He  went 
on  to  describe  a  method  for  choosing  the  order  using  lattice  filters.  Li  and 
Dickinson46  have  also  proposed  a  lattice  filter  approach  for  the  case  where  the 
noise  is  white. 

3.1.4  The  Maximum  Entropy  Method 

In  a  foundational  paper  on  information  theory,  Shannon47  described  the 
entropy  H  of  the  random  variable  x  as 

I 

H  =  -J^pilogpi  (3.2) 

1=1 

where  there  was  a  probability  p,-  that  the  random  variable  would  take  on  the 
value  x,.  In  a  paper  published  in  1957,  Jaynes48  proposed  that  the  the  under¬ 
lying  probability  density  functions  for  statistical  mechanics  should  be  modeled 
as  having  the  maximum  randomness,  or  maximum  entropy,  as  this  entailed  the 
least  number  of  assumptions.  This  concept  became  more  widely  known  in  the 
field  of  array  processing  as  the  maximum  entropy  method  (MEM)  due  to  Burg49 
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in  1967  when  he  applied  the  concept  to  spectral  estimation. 


A  compact  derivation  has  been  given  by  McDonough50  for  the  case  where 
the  noise  is  Gaussian.  For  a  uniformly  spaced  line  array  the  result  is 


B\0)  = 


Y  Vij  exp(j2nk(xi  —  xj)  cos  0 


1 


ET*VE 


(3.3) 


(3.4) 


where  the  variables  Vjy  are  derived  using  the  spatial  cross-correlation  matrix  3? 
and  the  requirement 


expf  j2nk(zi  —  zj)cosO 


— - * - L - dO 

£  Uiexpp2xfc(n  -  a;()  cos  (?) 

(3.5) 

eeT'  m 
et*ve 

(3.6) 

with  the  integral  extending  over  all  values  of  9  from  which  energy  is  arriving.  The 
symbols  ( T *)  denote  the  conjugate  transpose  operation  and  E  denotes  a  simple 
plane  wave  steering  vector.  As  noted  by  Johnson,25  Eq.  (3.3)  corresponds  to  the 
linear  predictor  or  autoregressive  solution;  thus,  all  results  from  the  previous 
discussion  of  the  AR  model  applies  here. 


3.1.5  Using  the  ARMA  Methods  to  Locate  Scatterers 

We  wish  to  consider  whether  these  methods  could  be  used  to  predict  the 
data  beyond  the  ends  of  a  line  array.  We  can  think  of  this  situation  as  being  one 
where  there  is  a  well-defined,  but  unknown,  scattering  function  which  describes 
the  strength  of  the  energy  reflected  from  the  source  or  sources  as  a  function  of 
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angle.  Our  finite-sized  line  array  has  recorded  a  subset  of  this  function,  over 
a  subset  of  angles.  The  questions  before  us  is  then  whether  we  can  use  this 
recording  to  predict  the  scattering  function  at  the  other,  unmeasured  angles. 

Such  a  prediction  would  seem  to  be  fundamentally  unsound  in  some  cases 
when  we  have  no  knowledge  of  the  underlying  scattering  function.  For  example, 
the  measured  scattering  patterns  from  spheres  shown  by  Pierce51  for  ka  >  1 
exhibit  a  complexity  near  the  nulls  that  would  be  hard  to  predict  from  a  recording 
of  only  a  small  subset  of  angles.  In  the  next  chapter,  Figure  4.1  shows  the 
diffraction  from  a  barrier;  note  that,  near  the  shadow  boundaries  at  0°  and  ISO0, 
there  are  fundamental  changes  in  the  pattern.  Likewise,  the  reflected  pattern  from 
plates  above  a  certain  size  exhibit  discontinuities  in  the  nature  of  the  pattern  at 
certain  angles,  as  is  shown  in  Figure  5.13.  No  empirical  method  for  determining 
the  correct  model  and  order  for  an  unknown  scatterer  is  known  to  this  author. 
It  would  also  follow  that  prediction  of  some  function  of  a  scattering  pattern  may 
also  be  difficult. 

A  fundamental  assumption  in  MEM  is  that  the  object  of  prediction  can  be 
modeled  as  having  a  Gaussian  distribution.  MEM,  in  effect,  expects  the  scattering 
pattern  as  a  function  of  angle  (or  its  Fourier  transform)  to  be  Gaussian  —  an 
unusual  if  not  impossible  case. 

Of  course,  if  we  are  sufficiently  far  from  the  source,  the  acoustical  wave 
is  merely  a  plane  wave  and  we  can  easily  extrapolate  it  to  larger  apertures  by 
shifting  the  phase  of  the  plane  wave.  However,  when  one  processes  the  data  from 
the  new,  simulated  receivers,  the  processing  function  must  also  be  phase  shifted 
by  the  same  amount;  thus,  the  result  is  essentially  the  same  as  multiplying  the 
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original  information  by  some  constant  and  then,  upon  normalization,  dividing  by 
the  same  constant. 

Such  considerations  seem  to  lead  to  the  conclusion  that  perhaps  the  ad¬ 
vantages  attributed  to  ARMA  by  Abdel-Aal  et  a/.43  are  actually  due  to  the  fact 
that  their  source  happened  to  match  the  ARMA  model. 

Another  problem  which  must  be  addressed  when  using  ARMA  (or  most 
other  high-resolution  methods)  is  their  tendency  to  produce  false  peaks  whenever 
they  model  the  source  field  as  an  AR  or  all-pole  process.  False  peaks  occurred 
when  Byrne  and  Fitzgerald52  used  MEM  with  an  array  whose  receivers  were 
spaced  less  than  a  half- wavelength  apart  although  Kaveh  and  Lippert53  suggested 
an  energy-taper  approach  that  eliminated  much  of  the  problem. 

3.2  Minimum  Energy  Methods 

In  the  1960’s  a  large  aperture  seismic  array  (LASA)  was  built  in  Montana 
to  measure  the  vector  velocity  of  propagating  seismic  waves.  The  purpose  of 
LASA  was  to  discriminate  distant  underground  nuclear  blast  tests  from  natural 
earth  tremors.  In  a  paper  published  in  1969,  Capon54  developed  a  high-resolution 
method  for  estimating  the  frequency-wavenumber  spectrum  received  by  this  ar¬ 
ray;  he  called  this  a  maximum  likelihood  (ML)  wavenumber  filter. 

3.2.1  Discussion  of  the  ME  Method 

This  method  is  an  extension  of  earlier  time  domain  work  by  Capon  et  al .55 
which  had  maximized  a  likelihood  function.  Johnson25  noted  that  the  method 
should  more  correctly  be  called  a  minimum  energy  (ME)  filter  since  it  minimizes 
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the  energy  in  a  beam  pattern  subject  to  a  constraint  but  does  not  maximize  a 
likelihood  function.  McDonough50  provided  an  elegant  method  for  deriving  this 
filter  that  is  based  upon  finding  a  wavenumber  filter  that  has  a  value  of  1  in 
the  direction  of  interest  6  but  which  otherwise  yields  the  minimum  energy.  His 
expression  for  the  ME  method  is 

BT>  -  («) 

where  E  is  a  simple  plane  wave  steering  vector  and  5R  is  the  spatial  cross¬ 
correlation  matrix  formed  from  the  zero-delay  correlation  between  all  possible 
receiver  pairs. 

This  beamformer  is  described  by  McDonough50  and  others  as  one  which 
reduces  its  side  lobes  the  most  in  those  directions  from  which  the  greatest  amounts 
of  energy  is  arriving;  which  allows  its  side  lobes  to  grow  in  those  directions  from 
which  lesser  amounts  of  energy  are  arriving;  and  which  passes  without  change 
plane  waves  from  the  direction  of  interest  0.  These  characteristics  have  led  to  the 
ME  method  being  described  as  a  wavenumber  filter.  Also,  because  it  does  not 
have  a  fixed  response  but  instead  varies  with  the  data,  it  must  also  be  classified 
as  a  non-linear  processor. 

Note  that  the  filter  itself  makes  no  distinction  between  signal  and  noise 
—  the  filter  attempts  to  suppress  energy,  regardless  of  its  source,  if  that  energy 
arrives  from  angles  other  than  the  one  of  interest.  When  operated  in  this  manner, 
the  ME  filter  is  said  to  be  based  upon  the  signal-plus-noise  cross-correlation 
matrix.  However,  when  it  is  desirable  and  possible  to  estimate  9?  from  signal-free 
inputs,  the  ME  filter  is  said  to  be  based  upon  the  noise-only  cross-correlation 
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matrix.  There  seems  to  be  no  other  standard  nomenclature  for  describing  these 
two  versions  of  the  ME  filter. 

3.2.2  Eigenvector  Analysis 

A  powerful  enhancement  to  the  ME  method  was  introduced56,57  when  it 
was  noted  that,  if  the  spatial  cross-correlation  matrix  3?  could  be  inverted,  it 
could  also  be  decomposed  into  a  set  of  orthogonal  eigenvectors.  Furthermore, 
the  magnitude  of  the  eigenvalue  associated  with  each  vector  corresponds  to  the 
energy  represented  by  the  vector.  When  the  number  of  energy  sources  is  limited 
to,  say,  iV5,  then  the  largest  Ns  eigenvectors  correspond  to  these  sources  and 
can  be  said  to  span  the  signal  subspace.  When  there  are  NR  receivers,  there  will 
then  be  NR  —  Ns  eigenvectors  with  relatively  small  eigenvalues  which  correspond 
to  noise  and  can  therefore  be  said  to  span  the  noise  subspace.  Since  eigenvectors 
are  mutually  orthogonal,  an  eigenvector  from  one  subspace  is  orthogonal  to  the 
other  subspace.  In  particular,  when  the  steering  vector  E  is  aimed  precisely  at 
the  bearing  of  an  actual  source,  it  will  reside  in  the  signal  subspace  and  so  it  too 
will  be  orthogonal  to  the  noise  subspace. 

This  characteristic  can  be  exploited  by  modifying  Eq.  (3.7)  so  that  31  is 
replaced  with  3Jre  which  represents  the  noise  subspace  and  is  constructed  using 
only  the  eigenvectors  and  eigenvalues  which  correspond  to  noise.  Then,  as  E 
approaches  a  true  source,  the  denominator  in  the  modified  Eq.  (3.7)  will  approach 
zero  and  the  beamformer  output  will  approach  infinity. 
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3.2.3  Linear  Predictors 

The  linear  predictor  method,  as  discussed  by  Johnson,25  is  based  upon 
the  prediction  of  a  selected,  single  element  of  the  array  by  a  linear,  weighted  sum 
of  the  other  elements  of  the  array: 


Rq  —  ~  amRm ■  (3.8) 

m^q 

The  true  value  of  Rq  is  known,  and  the  coefficients  am  are  found  by  minimizing 
the  mean-squared  error  subject  to  the  constraint  that  aq  =  1.  Johnson  notes  that 
this  is  a  constrained  optimization  problem  of  the  same  form  as  that  encountered 
for  the  minimum  energy  method,  and  also  corresponds  to  an  autoregressive  model 
of  the  signal.  The  solution  to  this  problem  is  then 


LPq(0)  = 


Uj*  3?-1  Uq 
{uj*  5R-1  E )2 


(3.9) 


where  Uq  is  a  column  vector  with  the  gth  element  equal  to  one  and  the  other 
elements  equal  to  zero.  The  name  of  this  function  includes  a  power  of  4  since  it, 
in  this  form,  is  proportional  to  the  square  of  the  power.  A  modified  version  is 
often  used  so  that  it  is  similar  to  other  beamformers  and  is  given  by 


£P,2(0)  = 


1 

l/T*  SR-1  £ 


(3.10) 


Typically,  the  value  for  the  the  center  receiver  or  one  of  the  receivers  at 
the  ends  of  a  line  array  are  chosen  for  prediction.  In  a  comprehensive  comparison 
between  the  LP,  ME,  and  Bartlett  beamformers,  DeGraff  and  Johnson58  found 
that  the  LP  method  generally  yielded  the  highest  resolution  although  it  needed 
large  amounts  of  averaging  and  could  also  produce  false  peaks  and  exhibit  in- 
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creased  bias  in  some  situations.  Furthermore,  Johnson20  notes  that  “a  criterion 
for  the  ‘proper'  choice  of  the  predictive  element  (q)  has  not  been  found  to  date.” 


3.2.4  Using  the  ME  Method  in  the  Absence  of  Noise 

Johnson25  has  expressed  the  inverse  of  the  cross-correlation  matrix  for  the 
case  of  white,  uncorrelated  noise  as 


ST1  = 


1 


I  - 


NR  <T2  +  <7“ 


EE 


T* 


(3.11) 


where  I  is  the  identity  matrix,  E  is  a  plane  wave  steering  vector,  and  the  power 
levels  of  the  signal  and  noise  are  given  by  a~  and  cr2,  respectively.  In  this  case, 
it  is  clear  that  letting  er2  — >  0  would  be  unacceptable.  The  eigenvector  analy¬ 
sis  suggests  why  this  occurs.  As  the  noise  goes  to  zero,  the  eigenvectors  of  3? 
corresponding  to  noise  become  identical,  their  eigenvalues  go  to  zero,  and  so  the 
inverse  of  no  longer  exists. 


Of  course,  even  a  simulation  will  suffer  noise,  even  if  only  digital  noise 
due  to  numerical  roundoff.  However,  it  is  not  clear  what  difficulties  may  be 
encountered  at  extremely  low  levels  of  noise;  in  these  cases,  some  of  the  techniques 
discussed  later  in  this  chapter,  for  accommodating  coherent  sources,  could  be 
employed. 


3.2.5  Using  the  ME  Method  to  Locate  Scatterers 

We  now  wish  to  explore  the  strengths  and  weaknesses  of  this  method  from 
the  point  of  view  of  one  wishing  to  use  a  linear  array  to  locate  scatterers. 
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3. 2. 5.1  The  ME  Wave  Model 

In  Section  2.2.2,  it  was  noted  that  the  Barlett  beamformer  could  fail  when 
its  model  of  the  energy  waves  (it  expects  plane  waves)  does  not  correspond  to 
the  actual  energy  wave.  The  ME  method  could  fail  for  a  similar  reason:  it  would 
interpret  a  complicated  wavefront  arriving  from  a  single  direction  as  simple  planar 
wavefronts  arriving  from  many  directions.  This  problem  has  been  addressed 
by  Seligson59  and  further  discussed  by  McDonough.60  The  problem  can  become 
so  severe  in  some  cases  that  the  ME  resolution  could  be  worse  than  that  of 
the  Bartlett  beamformer  or  even  degenerate  to  totally  useless  information.  The 
general  experience  has  been  that  the  problem  is  worst  when  there  are  a  relatively 
small  number  of  strong  sources  with  relatively  weak  background  noise. 

Seligson59  attempted  to  construct  a  variation  of  the  ME  method  which 
used  a  given  shape  for  the  wavefront.  He  found  that  the  processor  produced  worse 
resolution  than  the  Bartlett  beamformer  when  the  received  shape  differed  from 
the  expected  shape,  or  when  the  amplitudes  of  the  signals  differed  more  than  a 
certain  amount  between  receivers.  McDonough60  attempted  to  build  a  processor 
bat  was  more  tolerant  of  such  problems  by  extending  Seligson’s  approach  to 
consider  a  family  of  wave  shapes  consisting  of  limited  perturbations  of  the  given 
shape. 

The  central  ideal  of  the  ME  design  is  to  make  the  filter  transparent  to  all 
waves  from  the  direction  of  interest  d,  regardless  of  the  shape  of  the  wavefront. 
It  is  theoretically  possible  to  adapt  the  steering  vector  E  to  correspond  to  a 
wavefront  of  any  desired  shape;  however,  even  then,  two  problems  remain:  first, 
some  means  must  be  found  for  knowing  what  curvature  to  choose;  and  second, 
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the  filter  will  now  not  necessarily  be  the  minimum  energy  filter. 

3. 2. 5. 2  Number  of  Sources 

The  eigenvector  analysis  requires  that  the  number  of  sources  Ns  be  esti¬ 
mated  using  the  cross-correlation  matrix.  Although  several  methods  have  been 
proposed,61  it  remains  a  difficult  decision.  Johnson  and  DeGraaf62  found  that  the 
basic  ME  method  yielded  l\rS  peaks  no  matter  how  many  actual  sources  were  re¬ 
ally  present;  an  eigenvector  analysis,  on  the  other  hand,  was  much  more  tolerant. 
They  provided  one  example  where  three  peaks  were  produced,  correctly  locating 
the  three  actual  sources,  in  spite  of  having  set  Ns  —  6.  However,  they  conclude 
that  “reasonably  accurate  methods  of  determining  (the  number  of  sources)  from 
the  eigenvalues  of  5R  are  not  known  at  this  time.” 

Zou  and  Liu63  also  reported  false  peaks  using  the  LP  method  but  found 
that  improvements  could  be  made  by  using  the  average  of  two  LP  variations. 

The  eigenvector  analysis  can  accommodate  up  to  NR  —  1  uncorrelated 
sources  when  there  are  NR  receivers.  When  the  sources  are  correlated,  Bresler 
and  Macovski64  note  that  this  analysis  is  then  limited  to  \NR.  However,  they 
go  on  to  establish  a  theoretical  maximum  which  can  be  much  higher,  depending 
on  the  signal-to-noise  ratio  and  the  angular  distribution  of  sources.  Pillai  et  a/.40 
showed  that  NR(NR  —  1)  sources  could  be  resolved  if  the  receivers  were  spaced 
according  to  a  minimum-redundancy  pattern  called  a  Caratheodory  sequence. 

It  would  seem  that  an  eigenvector  analysis  would  have  difficulty  with 
an  extended  scatterer  which  acts  like  an  extremely  large  number  of  scattering 
sources.  However,  the  key  requirement  is  that  the  N *  eigenvectors  with  the 
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largest  eigenvalues  span  the  source  subspace  —  the  number  of  actual  or  equiva¬ 
lent  sources  does  not  itself  have  to  equal  Ns .  Thus,  the  results  depend  upon  the 
relative  geometry  of  the  array  and  the  scatterer.  A  simulation  by  Duckworth41 
was  indeed  able  to  provide  reasonable  results  for  a  particular  extended  scatterer, 
with  much  sharper  boundaries  than  a  conventional  beamformer.  The  general  con¬ 
ditions,  however,  under  which  extended  scatterers  can  be  accommodated  have  yet 
to  be  established. 

3. 2. 5. 3  Coherent  Sources 

When  two  or  more  of  the  sources  are  coherent,  we  can  find  neither  the 
inverse  nor  the  eigenvectors  of  the  cross-correlation  matrix  9ft.  This  will  result  in 
all  cases  when  the  sources  are  simply  scattering  energy  from  the  same  transmitter. 

A  method  called  spatial  smoothing  was  proposed  by  Shan  and  Ivailath65 
wherein  the  linear  array  was  treated  as  subarrays,  each  of  which  shared  all 

but  one  receiver  with  each  neighbor.  The  cross-correlation  matrices  from  each 
of  these  subarrays  were  then  averaged  to  form  a  single  cross-correlation  matrix 
in  which  the  coherence  between  sources  has  been  broken.  This  method  assumes, 
however,  that  the  full  array  be  such  that  it  can  be  subdivided  into  identical 
subarrays.  It  further  depends  upon  uncorrelated  noise  existing  at  the  receivers 
for  the  loss  of  coherence;  if  there  is  little  noise,  the  method  may  work  poorly. 
Worst  of  all,  the  effective  size  of  the  array  is  now  half  the  original  size  which 
works  against  the  high  resolution  we  are  attempting  to  achieve.  A  variation 
which  does  not  cause  as  great  a  reduction  in  the  effective  array  size  and  is  called 
the  modified  spatial  smoothing  (MSS)  method  was  described  by  Evans  et  al.66 
although  it  has  been  noted  that  the  method  may  fail  in  some  special  cases.67 
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A  similar  method  called  frequency  smoothing  was  proposed  by  Wang  and 
Ivaveh68  which  averages  the  3ft  obtained  at  each  frequency.  In  particular,  they 
combine  the  information  from  each  frequency  in  a  manner  they  describe  as  “coher¬ 
ent  addition”  so  as  to  avoid  the  threshold  effect  which  has  kept  similar  methods 
using  an  incoherent  addition  from  succeeding.  However,  as  noted  by  Williams  et 
al.,6‘  this  approach  requires  either  prior  knowledge  of  the  general  source  locations 
or  a  very  good  estimate. 

Duckworth41  found  that  sources  also  became  decorrelated  when  the  cross¬ 
correlation  matrices  from  several  sets  of  data  were  averaged  if  there  was  a  small 
amount  of  motion  between  the  transmitter,  source,  and  receiver  between  sets  of 
data. 


Several  other  methods  for  breaking  the  correlation  have  begun  to  be 
investigated69  including  a  method  by  Bresler  et  al.70  which  combines  each  group 
of  coherent  sources  into  one  equivalent  source. 
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Chapter  4 

Pattern-Match  Imaging 

It  was  claimed  in  Chapter  2  that  the  pattern-matching  approach  can  ac¬ 
commodate  events,  such  as  diffraction  and  reflection,  which  can  occur  during 
the  propagation  of  the  pressure  wave.  To  investigate  this  claim,  two  new  imag¬ 
ing  techniques,  called  the  pattern-match  method  and  the  mismatch  method,  are 
developed  in  this  chapter  using  the  concept  of  a  fictitious  source  introduced  in 
Section  2.4.  A  simulation  of  the  diffraction  of  waves  around  a  barrier  is  also 
implemented  and  is  used  to  examine  the  capabilities  of  the  new  methods.  Each 
simulation  will  also  be  repeated  using  the  well-known  holographic  reconstruction 
method  so  that  all  3  techniques  can  be  compared.  Additional  simulations  for  a 
selection  of  special  cases  will  be  included  to  round  out  the  investigation;  these 
cases  include  a  free-field  environment  (i.e.,  the  diffracting  barrier  is  eliminated), 
a  case  where  the  location  of  the  barrier,  used  in  the  reconstruction,  is  in  error; 
and  a  case  where  there  are  two  sources. 

4.1  Simulation  of  Edge  Diffraction 

The  performance  of  the  pattern-matching  method  depends  upon  the  de¬ 
gree  to  which  the  signals  RF  received  from  the  fictitious  source  agree  with  the 
true  physical  phenomena.  While  this  is  no  challenge  in  a  Iree-space  situation, 
it  may  become  a  limiting  factor  when  attempting  to  apply  the  method  to  cases 
where  the  pressure  wave  has  undergone  diffraction  and  reflection. 

When  the  received  signal  is  simulated,  as  it  will  be  in  this  chapter,  there 
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is  no  problem:  whatever  formula  we  use  to  simulate  the  received  signals  is  used 
to  represent  the  fictitious  source.  Nevertheless,  we  should  choose  a  non-trivial 
case  so  that  artifacts  and  limitations  of  the  method,  if  any,  might  be  observed. 
These  artifacts  may  include  bias  (the  source  location  is  in  error),  resolution  (the 
precision  with  which  the  source  has  been  located),  and  non-uniqueness  (false 
sources  are  reported). 

4.1.1  Defining  the  Wave  Disturbance 


The  scattering  of  pressure  waves  by  a  sphere  is  an  attractive  physical  phe¬ 
nomenon  for  research  since  a  entire  series  of  scattering  patterns  can  be  generated 
merely  by  adjusting  the  term  ka  where  k  is  the  wavenumber  and  a  is  the  radius 
of  the  sphere.  Pierce71  writes  the  scattered  component  as 


Psc  — 


exp(j  k  r) 
r 


(4.1) 


where  B  is  the  peak  amplitude  of  the  incident  planp  wavp,  r  is  the  distance  from 
the  sphere  to  a  point  in  the  scattered  field,  and  0  is  the  angle  between  r  and  the 
direction  of  the  incident  wave.  Unfortunately,  this  simple  formula  holds  only  for 
ka  <C  1  which  translates  to  a  ~  0.016A  in  this  study  where  the  wavelength  is  one. 
This  radius  enters  into  Eq.  (4.1)  in  the  term  that  calculates  the  volume  of  the 
sphere  and  results  in  the  scattered  component  having  an  amplitude  five  to  six 
orders  of  magnitude  lower  than  the  non-scattered  component.  It  was  estimated 
that  when  considering  the  total  field,  equal  to  the  sum  of  these  two  components, 
the  scattered  component  amplitude  would  be  too  small  an  effect  to  demonstrate 
the  abilities  of  pattern  matching  over  holography. 


A  more  pronounced  effect  would  be  generated  by  the  edge  of  a  barrier. 
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Pierce' “  has  expressed  the  diffracted  component  of  the  field  when  the  distance  of 
both  the  source  and  field  point  from  the  edge  is  large  compared  to  a  wavelength 
as 


exP  (j  k  (r3  +  ry)  1  exp  (jf) 
Pd~b  rs  +  rf  72 


E 


sin(^  7r) 


(4.2) 


Ar:(r  ■  A/„  (9- 


(4.3) 


+.-  A  —  cos(i/  7 r)  •  cos  (u  ( 
where  the  summation  consists  of  two  terms.  The  angle  comparison  <p discussed 
further  in  the  next  section,  is  given  by 

(  <f>R  +  <f>s  +  37 r,  if  x  is  plus; 

4>x  =  s 

(  < f>R  —  <})S ,  if  x  is  minus. 

The  source  is  located  a  distance  ra  from  the  edge  of  the  barrier  at  angle  <f)S  and 
transmits  waves  of  amplitude  S.  The  field  point  where  we  wish  to  evaluate  the 
diffracted  field  is  a  distance  rj  from  the  edge  at  angle  <f>R.  The  wedge  index  v  is 
^  for  a  thin  barrier.  We  also  have 


r  = 


\ 


2  77  rs 


MM)  = 


A  (77  +  rs)  ’ 

COs(l/7r)  —  COs(l/  <$>) 
U\J  1  —  cos(j/  7 r)  •  cos(u  <j>) 

— 2cos(j<^>)  for  v  = 


(4.4) 

(4.5) 

(4.6) 


For  simplicity,  we  set  A  =  1.  The  term  Ap  is  the  diffraction  integral  which  can 
be  written  in  terms  of  the  Fresnel  integrals 


C(fi)  =  j  cos  (f  t2)  dt, 
S(f 1)  =  J  sin  (?f2)  dt. 


(4.7) 
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as 


Ad(v)  =  ^-~exp  (-;  §  V2)  (sign (ft)  -  (1  -  j)(c(n)  +  jS(fi))')  (4.8) 


where 


SIGN(/x)  = 


1,  fi>  0; 


(4.9) 


—  1,  /LX  <  0. 

Some  routines  for  the  evaluation  of  the  Fresnel  integrals  use  the  absolute  value 
of  their  argument,  thereby  failing  to  include  the  odd  symmetry.  In  these  cases, 
we  can  move  the  SIGN  operation  to  correct  for  this  so  that  Eq.  (4.S)  becomes 

Ap(/x)  =  exp  (-j  f  n2)  sign (n)  (l  -  (1  -i)(C(|/i|)+j5(|/i()))  (4.10) 

where  | . . .  j  denotes  the  absolute  value  operation. 


4.1.2  Angle  Definitions 

It  is  convenient  to  define  the  angles  which  describe  the  location  of  the 
source  (<j>s)  and  field  point  (<f>R)  in  the  manner  most  useful  for  the  situation. 
Unfortunately,  this  has  led  to  three  different  reference  points  for  these  angles  and 
so  a  note  on  these  definitions  is  appropriate.  In  all  cases,  the  angles  increase  in  a 
counter-clockwise  direction  —  these  definitions  differ  only  in  the  location  of  the 
0C  point. 


On  polar  plots,  the  angular  origin  will  be  on  the  positive  Z-axis,  as  is 
shown  in  Figure  4.1,  and  will  be  expressed  in  degrees.  In  this  figure,  two  angles 
are  shown  as  examples:  <ps ,  the  angle  of  the  source;  and  <f>R ,  the  angle  of  a  receiver 
or  field  point.  Note  that  both  of  these  are  examples  of  negative  angles. 


The  radial  distance  of  the  source  on  this  and  subsequent  figures  is  not 


source 


Figure  4.1.  Regions  of  the  field  scattered  by  a  barrier,  numbered  according  to  the  number  of 
components  used.  The  location  of  the  regions  will  change  if  the  location  of  the  source  changes. 
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drawn  to  scale  —  in  fact,  on  subsequent  plots,  radial  distances  will  instead  rep¬ 
resent  the  signal  level  in  dB. 

The  sketch  of  the  barrier  is  also  not  drawn  to  scale;  it  is  included  at  the 
+90°  location  only  as  a  reminder  of  the  angular  location  of  the  barrier. 

In  Figure  4.1,  the  edge  of  the  barrier  is  shown  as  being  located  at  the 
origin  of  the  coordinate  system.  This  convention  will  be  used  on  all  subsequent 
polar  plots.  However,  in  the  simulations,  the  edge  of  the  barrier  is  never  at  the 
origin  but  is  instead  at  some  other  given  location.  This  requires  a  coordinate 
transformation,  which  is  made  automatically  in  the  subroutine  which  calculates 
the  diffraction. 

The  second  angular  reference  point  is  found  in  the  fortran  routines  where 
it  was  more  convenient  to  put  the  angular  origin  on  the  positive  X-axis  so  that 
the  inverse  tangent  routine  could  directly  use  the  X  and  Z  coordinates  of  a  point 
in  the  calculation  of  the  angle.  The  third  angular  reference  point  is  found  in 
Pierce71  in  his  development  of  the  diffraction  equations  where  he  sets  the  angular 
origin  at  the  shadow  side  of  the  barrier,  which  would  be  the  negative  X  axis  in 
this  study. 

The  angles  used  on  the  plots  will  agree  with  those  used  in  equations  in 
this  study  so  that  the  reader  need  not  consider  these  other  definitions  except 
when  examining  the  computer  code  or  other  texts.  The  necessary  conversion  is 
embedded  in  the  offset  of  added  to  each  angle  in  Eq.  (4.3). 
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4.1.3  The  Total  Field 

The  total  field  consists  of  the  diffracted  component  (pj),  discussed  previ¬ 
ously,  plus,  in  some  regions,  a  directly  incident  component  (p,)  and  a  reflected 
component  ( pT ).  These  regions  can  be  uniquely  numbered  according  to  the  num¬ 
ber  of  components  which  are  needed  since  all  regions  include  p^  while  none  include 
pr  without  pi . 

Figure  4.1  shows  the  regions  which  result  for  a  particular  source  location. 
The  line  dividing  Region  1  from  Region  2  is  called  the  incident  shadow  boundary, 
while  the  line  dividing  Region  2  from  Region  3  is  called  the  reflection  shadow 
boundary. 

Table  4.1  shows  the  logic  which  is  used  in  the  subroutine  to  determine  the 
region  number  based  upon  the  definition  of  $  used  in  the  plots  and  equations. 


Table  4.1.  Rules  for  determining  the  region  number. 
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The  incident  component,  when  needed,  is  given  by 


Pi  =  S 


exp (j  k  r) 


(4.11) 


where  r  is  the  direct  distance  between  the  source  and  the  field  point.  The  reflected 


component,  when  needed,  is  given  by 


pr=S 


exp  J 


1  k  (rxa  + 


rx a  ■  rxf 


(4.12) 


where  rx3  and  rxf  are  the  distances  between  the  point  of  reflection  and  the  source 
and  field  point,  respectively.  Given  the  source  location  (Xs ,  Zs),  the  field  point 
location  (XF,  ZF),  and  the  barrier’s  coordinate  ZB  (it  has  no  X  coordinate  as  it 
is  parallel  to  the  X-axis),  the  point  of  reflection  (X,  ZB)  is  given  by 


XF  •  \ZB  -  Zs |  +  Xs  •  \ZB  -  ZF\ 
\ZB  -ZS\-\ZB  -ZF\ 


(4.13) 


The  sign  of  the  complex  exponents  in  Eq.  (4.11)  and  (4.12)  are  positive  in 
order  to  agree  with  Eq.  (4.2).  However,  the  complex  conjugate,  of  the  sum  of  all 
components,  must  be  taken  so  that  the  phase  agrees  with  the  definition  of  phase 
used  in  this  study. 


4.1.4  Shadow  Boundary  Approximation 


When  the  receiver  is  near  the  shadow  boundary,  Pierce'2  has  developed 
an  approximate  solution  which  would  replace  Eq.  (4.2)  with 


Pd  =  s 


exp (jfc(r,  +  r/))  exp(jf) 


r3  +  rf 


AD(rr)- 


(4.14) 


This  approximation  could  cut  the  number  of  evaluations  of  the  diffraction  integral 
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Ad  in  half  and  so  will  be  considered  in  the  next  section. 

4.1.5  Testing  the  Barrier  Simulation 

The  absolute  value  of  the  Fresnel  integrals  were  calculated  in  a  FORTRAN 
subroutine  using  a  power  series.  Results  for  a  selection  of  arguments  from  0.0 
to  5.0  were  compared  to  the  tables  in  Abramowitz  and  Stegun.'3  The  values 
differed  by  at  most  ±1  in  the  last  digit;  that  is,  the  computed  results  seemed 
to  be  accurate  to  ±10-7.  Single-precision  arithmetic  (i.e.,  32  bits  were  used  to 
represent  a  number)  caused  2-5  times  as  many  differences  between  the  subroutine 
and  the  tables;  therefore,  double- precision  arithmetic  (using  64  bits  to  represent 
a  number)  is  used  in  the  subroutine. 

Another  fortran  subroutine  then  makes  use  of  the  Fresnel  integrals  to 
evaluate  the  diffraction  integral  and  the  pressure  at  given  points.  This  subroutine 
was  tested  by  evaluating  the  diffracted  wave  for  the  case  of  a  source  at  ( X ,  Z)  = 
(0, 100A)  and  a  barrier  beginning  at  (0, 10A)  and  extending  towards  X  =  — oo. 
The  magnitude  and  phase  of  the  field,  at  points  located  every  quarter  degree  on  a 
circle  of  radius  10A  around  the  edge  of  the  barrier,  is  shown  in  Figures  4.2  and  4.3, 
respectively.  The  phase  shown  has  been  “linearized”  by  adding  or  subtracting 
integer  multiples  of  360°  to  each  value  so  as  to  minimize  the  difference  between 
the  value  and  its  neighboring  values.  These  results  look  reasonable  and  seem 
compatible  with  a  similar  case  described  by  Kendig  and  Hayek.74 

The  same  case  using  the  shadow  boundary  approximation  is  shown  in  Fig¬ 
ure  4.4.  Given  the  assumptions  which  were  used,  the  results  are  remarkably  good. 
Nevertheless,  an  increase  in  the  relative  size  of  the  ripples  at  negative  angles,  and 


source 
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was  evaluated  at  a  radius  of  10A  from  the  edge. 


source 


evaluated  at  a  radius  of  10\  from  the  edge. 


source 
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Figure  4.4.  Magnitude,  using  the  shadow--boundary  approximation,  of  the  pressure  from  a  source  100A 
away  from  the  edge  of  a  barrier.  The  field  was  evaluated  at  a  radius  of  10A  from  the  edge. 
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a  discontinuity  at  0  =  0°,  would  make  this  approximation  unacceptable  for  some 
situations.  Testing  has  also  indicated  that  the  time  saved  not  significant  for 
this  study,  and  so  this  approximation  was  not  used  further. 

4.2  Mismatch  Imaging 

During  the  investigations  of  this  chapter,  it  was  noted  that  the  pattern- 
match  method  seemed  to  yield  a  poorer  image  in  some  situations  than  was  ex¬ 
pected.  That  is,  the  value  of  z\  ■  z%,  where  *  denotes  the  complex  conjugate, 
decreases  more  slowly  than  was  desired  as  the  complex  numbers  z\  and  zo  began 
to  differ. 


An  approach  to  improving  this  situation  is  suggested  by  the  minimum 
energy  (ME)  method  discussed  in  the  previous  chapter.  In  that  method,  the 
“match”  is  provided  by 


B2(0) 


1 

ER~lE 


(4.15) 


where  E  is  a  vector  that  steers  the  reception  to  angle  0  and  as  such  somewhat 
plays  the  role  of  the  fictitious  source  in  the  pattern-match  method.  The  basic 
form  of  Eq.  (4.15)  suggests  that  the  denominator  is  calculating  the  amount  of 
mismatch,  rather  than  match,  between  the  fictitious  source  of  plane  waves  and 
the  patterns  actually  received.  Thus,  when  the  fictitious  source  approaches  an 
actual  source,  the  amount  of  mismatch  approaches  zero  and  the  sharp  peaks 
typical  of  the  ME  method  result. 


What  we  desire,  therefore,  is  a  comparison  of  two  complex  values  which 
yields  a  result  proportional  to  the  degree  to  which  they  differ.  In  the  complex 
plane,  every  complex  value  is  represented  by  a  point,  and  so  the  concept  of  the 
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distance  between  any  two  points  z\  and  zo  could  be  defined  as 

zd  =  \f(a\  -  ao)~  +  (61  -  bo)2  (4.16) 

where  z\  =  (cti  +  jbi)  and  zo  =  («2  + j^)-  Actually,  we  desired  z~t  so  that  a 
squared  quantity  results  as  with  z\  ■  z%  or  the  correlation  between  two  arrays 
of  numbers.  This  quantity  will  be  used  to  implement  what  shall  be  called  the 
mismatch  method.  Any  comparison  method  using  z\  ■  zX  shall  continue  to  be 
called  the  pattern-match  method.  In  subsequent  tests  in  this  chapter,  the  results 
of  both  match  methods  will  be  compared  to  that  of  holography. 

In  plots,  it  is  conventional  for  higher  values  to  represent  stronger  peaks. 
We  need,  therefore,  to  either  use  \/z\  or  — z The  former  would  exhibit  the  non¬ 
linear  response  typical  of  the  ME  method  and  would  actually  be  quite  effective 
for  finding  point  sources.  However,  we  choose  to  use  the  latter  because  it  retains 
its  linearity  and  can  therefore  be  more  directly  compared  to  the  results  from  the 
pattern-match  method  and  holography. 

Using  the  pattern-match  method,  the  lower  bound  of  zero  represents  the 
worst  possible  match  while  the  upper  bound,  representing  the  best  match,  can 
be  easily  normalized  to  a  peak  value  of  one.  We  wish  to  map  the  new  values 
from  the  mismatch  method  to  this  same  domain.  The  upper  bound  of  —  zj  is 
zero,  representing  the  best  possible  match  (or  worst  possible  mismatch),  while 
the  lower  bound  is  some  negative  number  corresponding  to  the  worst  match  (or 
greatest  mismatch).  While  it  is  clear  that  these  values  should  be  adjusted  so  that 
the  best  match  is  always  mapped  to  one,  it  is  not  so  clear  which  (negative)  value 
ought  to  be  mapped  to  zero.  After  some  experimentation,  the  most  useful  results 
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were  obtained  by  mapping  the  average  (negative)  value  to  zero,  and  setting  values 
less  than  the  average  (representing  very  poor  matches)  equal  to  zero. 

This  transformation  is  implemented  by  first  calculating  all  values  of  z;j, 
their  average  zj,  and  the  minimum  value  corresponding  to  the  best  match  zj. 
Then,  the  result  of  the  new  method  $  which  will  be  used  for  plotting  and  com¬ 
parison  is  given  by 

T(x,2)  =  max  f  0,  - d  (4-17) 

V  zd  -  zd  J 

where  the  max  operation  chooses  the  larger  of  its  two  operands. 

4.3  Normalization 

In  early  tests,  the  pattern  match  method  failed  to  produce  a  maxima  in 
some  cases  when  the  fictitious  source  passed  through  the  location  of  the  true 
source.  To  investigate,  one  particular  case  was  examined  in  depth.  In  this  case, 
a  barrier  beginning  at  (0, 10A)  and  extending  towards  X  =  —00  diffracts  the 
pressure  waves  from  a  point  source  at  (— 5.25A,  15A).  An  array  of  64  receivers 
spaced  0.5A  apart  is  centered  on  the  X-axis.  The  pattern  received  at  this  array, 
and  the  pattern  from  a  fictitious  source  at  (— 4.75A,  15A),  are  shown  in  Figure 
4.5.  When  the  fictitious  source  was  at  the  exact  location  of  the  source,  its  pattern 
matched  exactly  the  received  pattern  (the  solid  line  in  Figure  4.5). 

When  the  fictitious  source  was  at  the  true  source,  the  match  value  was 
0.0162.  However,  when  the  fictitious  source  moved  to  the  location  shown  in 
Figure  4.5,  the  match  value  increased  to  0.0183.  This  was  the  phenomenon  which 
was  causing  the  pattern-match  method  to  fail  to  produce  the  expected  maxima. 

What  has  happened  is  that  the  pattern  from  the  fictitious  source,  as  it 
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moved  towards  the  edge,  increased  enough  (more  energy  was  getting  around  the 
barrier)  to  more  than  compensate  for  the  misalignment  of  the  patterns.  This 
phenomenon,  it  turns  out,  is  often  observed  in  signal  processing,  and  can  be 
elegantly  suppressed  by  normalizing  the  power  of  each  pattern  before  comparison: 

*  -  m  (4'1S) 

Thus,  the  absolute  and  relative  amplitudes  are  changed,  as  shown  in  Figure  4.6. 

The  effects  of  this  normalization  can  be  dramatic.  The  results  shown  in 
Figure  4.29  later  in  this  chapter  had  originally  appeared  as  shown  in  Figure  4.7 
when  the  pattern-match  method  was  used  without  normalization  (note  that  these 
figures  differ  in  their  rotation  by  180°).  The  rapid  rise  near  Z  —  0  in  Figure  4.7 
is  due  to  the  pi  component  whose  f  term  grows  to  infinity  as  the  fictitious  source 
approaches  the  location  of  a  receiver. 

The  large  positive  peaks  in  Figure  4.7  become  large  negative  peaks  when 
the  mismatch  method  is  used.  Thus,  the  results  shown  in  Figure  4.33  later  in 
this  chapter  had  originally  appeared  as  shown  in  Figure  4.8  without  normaliza¬ 
tion.  The  elevated  area  or  plateau  is  due  to  the  large  negative  peaks  causing  an 
inappropriate  shift  in  the  value  of  z\  used  in  Eq.  (4.17). 

4.4  Locating  a  Source  in  Free  Space 

A  simulation  was  prepared  of  a  single  harmonic  point  source  whose  sig¬ 
nals  were  recorded  at  an  array  of  64  receivers.  The  point  source  was  located 
at  (X,Z)  =  (0.25A,15A).  This  slightly  off-center  point  was  chosen  so  that  the 
source  would  fall  upon  one  of  the  reconstruction  locations;  this  prevents  some  of 
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ire  4.7.  Reconstruction  using  the  pattern-match  method,  without  normalization,  of  a  point  source 
-5.25A.15A)  with  diffraction  by  a  barrier  extending  from  (0, 10A)  to  (-<*>, 10A). 


so 
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Figure  4.8.  Reconstruction  using  the  mismatch  method,  without  normalization,  of  a  point  source  at 
(-5.25A.15A)  with  diffraction  by  a  barrier  extending  from  (0,10A)  to  (-<*>, 10A). 


SI 

the  minor  anomalies  which  can  occur  when  the  source  falls  between  reconstruc¬ 
tion  locations.  The  receivers  were  located  on  the  A  axis  at  intervals  of  AA  with 
the  array  centered  so  that  its  ends  were  at  ±15.7oA. 

In  this  section,  and  all  subsequent  sections  in  this  chapter,  the  selection 
of  the  Z  coordinate  of  points  in  the  reconstructed  field  is  arbitrary.  Thus,  we  are 
free  to  select  values  which  span  the  region  of  interest,  i.e.,  from  near  the  receiver 
array  (Z  =  1A)  to  beyond  the  source  location  (typicallv,  Z  =  15A  or  Z  —  20A). 
The  reconstructed  field  was  evaluated  at  intervals  of  1  wavelength  (i.e..  Z  =  1A, 
Z  =  2A,  . .  .  )  so  that  it  was  unlikely  that  any  feature  would  be  missed. 

The  choice  of  the  X  coordinate  of  points  in  the  holographic  reconstructions 
is,  unfortunately,  not  arbitrary  —  the  reconstructed  points  must  lie  at  the  same 
A’  coordinates  as  the  receivers  due  to  the  demands  of  the  FFT.  Although  this 
restriction  does  not  apply  to  the  pattern-match  or  mismatch  methods,  it  was 
adhered  to  in  all  of  the  reconstructions  produced  in  this  chapter  so  that  the 
results  could  be  more  easily  compared. 

In  the  parlance  of  holography,  the  field  points  are  the  locations  at  which 
the  pressure  field  was  “reconstructed.”  In  the  parlance  of  pattern  matching,  the 
fictitious  source  was  placed  at  these  locations.  In  this  stud)-,  the  word  “recon¬ 
structed”  will  henceforth  be  used  to  refer  to  the  estimate  of  the  source  location(s) 
regardless  of  the  method  by  which  that  estimate  is  made. 

The  holography  method  was  implemented  in  a  fortran  program,  while 
both  pattern-matching  methods  were  implemented  in  a  second  fortran  program. 
Input  parameters  in  these  programs  selects  whether  or  not  the  fictitious  source 


is  obstructed  by  a  barrier:  whether  or  not  the  shadow-boundary  approximation 
is  used:  and  which  match  method  is  used. 

Figure  4.9  exhibits  the  expected  curve  when  the  holography  method 

reconstructed  the  held  at  the  source  (i.e..  at  Z  =  15A). 

In  Figure  4.10,  equal-level  or  contour  lines  through  the  field  have  been 
generated,  while  in  Figure  4.11,  a  hidden-line  surface  plot  has  been  used  to  present 
the  results.  The  contour  lines  use  a  linear  interpolation  between  adjacent  .Y  values 
and  adjacent  Z  values  to  improve  the  quality  of  the  plot.  While  this  does  not 
distort  the  data,  it  may  produce  artificially  sharp  bends  in  some  lines.  Also,  an 
occasional  spike  in  the  contour  lines  may  result  when  the  algorithm  attempts  to 
determine  a  contour  through  a  group  of  points  which  have  nearly  the  same  value, 
such  as  near  ( —  1 A ,  5 A )  in  Figure  4.10. 

The  reconstruction  using  the  pattern-match  method,  shown  in  Figure  4.12, 
is  virtually  identical  with  Figure  4.11  —  only  a  small  reduction  near  Z  —  0  can 
be  seen.  The  reconstruction  using  the  mismatch  method  was  somewhat  difficult 
to  display  as  it  consisted  of  many  vertical  spikes.  Thus,  in  Figure  4.13,  only  a 
portion  of  the  information  around  the  source  is  shown,  and  the  values  which  fell 
below  15%  of  the  maximum  have  been  erased  from  the  plot.  In  Figure  4.14,  this 
same  portion  of  the  information  is  shown  as  a  contour  plot. 

One  (of  many)  interpretations  of  a  reconstruction  is  that  it  is  estimating, 
at  each  reconstructed  location,  the  “probability”  that  there  is  a  source  at  that 
location.  Ideally,  by  using  better  reconstruction  methods,  we  can  “push  down” 
the  probability  at  as  many  reconstruction  locations  as  is  possible;  this  helps 
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eliminate  these  locations  horn  consideration,  leaving  (hopefully)  only  one  location 
where  the  source  could  possibly  be  located. 

These  figures,  and  many  others  not  shown  here,  lead  to  the  observation 
that  the  envelope  of  the  mismatch  spikes  is  quite  similar  to  the  results  provided 
by  the  holography  and  pattern-match  methods.  Thus,  the  mismatch  method, 
at  the  location  of  its  spikes,  is  making  the  same  estimate  as  the  other  methods 
—  no  better,  no  worse.  On  the  other  hand,  at  the  nulls  between  the  spikes, 
the  mismatch  method  has  pushed  down  the  probability,  thus  eliminating  these 
regions  from  consideration.  In  particular,  measurements  of  Figure  4.14  indicate 
that  the  peaks  around  the  true  source  location  form  partial  rings  separate  by 
about  one  wavelength.  This  seems  reasonable  since  the  transmitter  is  producing 
sinusoids  which  exactly  repeat  every  wavelength. 

4.5  Locating  a  Source  Near  a  Reflector 

A  barrier  is  now  introduced  which  extends  from  (X,  Z )  =  (0, 10A)  to 
(—oo,  10A).  Otherwise,  the  simulation  is  identical  with  that  of  the  previous  sec¬ 
tion.  The  barrier’s  height,  in  the  Y  direction,  is  irrelevant  since  the  simulation  is 
being  conducted  in  only  two  dimensions.  The  field  generated  by  a  source  in  this 
location  is  shown  in  Figure  4.15. 

The  holographic  reconstruction,  shown  in  Figure  4.16,  exhibits  a  false 
source,  as  would  be  expected,  due  to  the  reflection  of  the  true  source  by  the 
barrier. 

The  pattern-match  reconstruction,  shown  in  Figure  4.17,  does  not  exhibit 
a  false  source  but  does  show  some  additional  peaks.  These  additional  peaks  can 


Holographic  recons Iruction  of  a  point  source  at  (-5.25A,7.£ 


Pattern-match  reconstruction  of  a  point  source  at  (-5.25A.7.5A)  with  a  barrier  extending 
to  (— °°,10A);  contour  lines  are  at  10%,  20%,  90%  of  the  maximum. 
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be  seen  more  clearly  in  Figure  4. IS  where  only  the  information  near  the  true 
source  has  been  retained. 

The  mismatch  method  reconstruction,  shown  in  Figure  4.19,  has  had  val¬ 
ues  below  10%  of  the  maximum  removed  for  clarity.  A  contour  plot  of  this 
reconstruction,  shown  in  Figure  4.20,  exhibits  partial  rings  of  peaks  around  the 
true  source,  separated  by  about  one  wavelength,  similar  to  that  shown  in  Figure 
4.14  for  the  free-space  case. 

In  addition,  the  mismatch  method  exhibits  some  low-level  matches  along 
lines  behind  the  barrier  with  the  lines  separated  by  about  a  wavelength.  There 
may  be  several  reasons  for  these  matches  behind  the  barrier.  First,  in  this  region, 
one  of  the  basic  assumptions  of  the  diffraction  equation  has  been  partially  violated 
since  some  of  the  points  are  not  many  wavelengths  from  the  edge.  Second,  as 
shown  by  Figure  4.6,  the  signal  at  the  receiving  array  from  points  behind  the 
barrier  look  very  much  like  a  sine  wave  and  so  may  produce  a  low  level  match 
with  almost  anything.  Such  a  phenomenon  would  somewhat  repeat  for  every 
additional  wavelength  that  the  fictitious  source  is  moved  away  from  the  edge  of 
the  barrier. 

In  many  imaging  situations,  the  sources  are  identified  as  being  at  those 
locations  where  the  image  level  is  above  some  threshold.  If  a  threshold  were  set 
at,  say,  50%  of  the  maximum,  the  mismatch  method  would  provide  the  most 
precise  location  of  the  source.  At  a  lower  threshold,  all  methods  exhibit  some 
problems:  both  matching  methods  show  additional  peaks;  and  the  holography 
method  shows  a  false  source  due  to  reflection. 
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Figure  4.18.  Pattern-match  reconstruction  of  a  point  source  at  (-5.t25A.7-5A)  with  a  barrier  extending 
from  (0.10A)  to  ( — 10A) . 


re  4.19.  Mismatch  reconstruction  of  a  point  source  at  (-5.25A.7.5A)  with  a  barrier  extending  from 
IA)  to  f-oo.lOA).  Values  below  10%  of  the  maximum  have  been  removed  for  clarity. 


Lch  reconslruclion  of  a  point  source  at  (-5.25A.7.5A)  with  a  barrier  extending  from 
contour  lines  are  at  10%,  20%,  ....  90%  of  the  maximum. 
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4.6  Locating  a  Source  Behind  a  Barrier 

The  situation  investigated  in  the  previous  section  is  now  repeated  except 
that  the  source  is  moved  behind  the  barrier  so  that  the  transmitted  energy  is,  to 
various  degrees,  blocked  by  the  barrier.  The  source  will  located  at  one  of  three  X 
locations  along  the  line  at  Z  =  15A:  X  =  0.25A;  X  =  —  2.25A;  and  X  =  —  5.25A. 
The  field  for  the  first  and  last  of  these  three  locations  is  shown  in  Figures  4.21 
and  4.22,  respectively. 

The  reconstructions  provided  by  the  holography  method  are  shown  in 
Figures  4.23-4.25.  A  comparison  with  the  free  space  holographic  reconstruction 
shows  that,  as  the  source  becomes  hidden  behind  the  barrier,  the  reconstructed 
peak  moves  closer  to  the  edge  of  the  barrier.  When  the  source  is  almost  com¬ 
pletely  hidden  (Figure  4.25),  the  peak  at  the  barrier's  edge  seems  to  be  the  only 
disturbance  —  there  is  little  to  even  suggest  that  there  is  a  source  behind  the 
barrier. 


The  pattern-match  reconstruction  for  a  source  at  (0.25A,  15A)  is  shown  in 
Figure  4.26  and  seems  to  be  a  little  poorer  than  the  corresponding  holographic 
reconstruction  (Figure  4.23)  because  the  main  peak  seems  to  have  become  slightly 
spread.  In  addition,  Figure  4.27  shows  that  there  is  a  slightly  higher  match 
in  the  “hidden  region”  which  falls  between  the  barrier  and  what  will  be  called 
the  boundary  line.  This  boundary  line  begins  at  the  edge  of  the  barrier  and 
runs  through  the  true  source  location  onto  infinity.  As  the  source  is  moved  to 
locations  further  behind  the  barrier,  the  main  change  in  the  reconstruction,  shown 
in  Figures  4.2S  and  4.29,  is  that  the  match  in  the  hidden  region  rises. 


source 


Figure  4.21.  Magnitude  of  the  pressure  from  a  source  at  (0.25A.15A)  with  a  barrier  extending  from 
(0.10A)  to  (— *»,10A).  The  field  was  evaluated  at  a  radius  of  10A  from  the  edge. 
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Figure  4.26.  Patlern— match  reconstruction  of  a  point  source  at  (0.25A.15A)  with  a  barrier  extending 
from  (0,10A)  to  (— °°,10A);  contour  lines  are  at  10%,  20%,  ...,  90%  of  the  maximum. 
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Figure  4.27.  Pattern— match  reconstruction  of  a  point  source  at  (0.25A.15A)  with  a  harrier  extending 
from  (0,10A)  to  (-°°,10A). 
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Figure  4.28.  Pattern-match  reconstruction  of  a  point  source  at  (-2.25A.15X)  with  a  barrier  extending 
from  (0.10A)  to  (  °°,  10X) . 
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Figure  4.29.  Pattern-match  reconstruction  of  a  point  source  at  (-5.25A.15A)  with  a  barrier  extending 
from  (0,10A)  to  (-°°,10A). 
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The  mismatch  reconstruction  for  a  source  at  (0.25A,  15A),  shown  in  Figure 
4.30.  exhibits  peaks  which,  as  in  earlier  cases,  have  an  envelope  similar  to  the 
pattern-match  reconstruction  and  a  spacing  roughly  equal  to  a  wavelength.  As 
seen  in  Figure  4.31,  the  peak  at  the  location  of  the  true  source  is  the  largest, 
but  only  by  a  small  amount.  Likewise,  as  the  source  is  moved  further  behind 
the  barrier,  Figures  4.32  and  4.33  exhibit  peaks  falling  within  the  envelope  of  the 
corresponding  pattern-match  results. 

These  cases  suggest  that  some  general  rules  might  be  stated  which  roughly 
describe  how  each  imaging  method  responds  as  the  source  becomes  hidden  fur¬ 
ther  behind  the  barrier:  (1)  holographic  imaging  moves  it’s  reconstructed  peak 
towards  the  edge  of  the  barrier;  (2)  pattern-match  imaging  indicates  only  that 
the  source  is  somewhere  in  the  hidden  region,  with  a  decreasing  likelihood  that 
the  source  is  near  the  boundary  line  as  the  true  source  moves  further  behind  the 
barrier;  and  (3)  mismatch  imaging  produces  peaks  within  the  envelope  of  the 
pattern-match  results,  arranged  in  groups  which  are  multiples  of  a  wavelength 
from  the  barrier  edge. 

4.7  Environmental  Errors 

The  two  pattern-matching  reconstruction  methods  use  knowledge  of  the 
barrier  in  their  analysis  of  the  received  signals.  It  would  not  be  unusual  for  there 
to  be  small  errors  in  this  knowledge;  therefore,  the  case  of  a  source  behind  the 
edge  of  a  barrier  was  repeated  with  an  error  deliberately  introduced.  The  results 
for  the  pattern-match  and  mismatch  methods  are  shown  in  Figures  4.34  and  4.35, 
respectively.  In  these  simulations,  the  received  signals  were  calculated  using  the 
true  location  of  the  barrier,  at  Z  —  10A,  while  the  reconstruction  used  a  false 


[ismatch  reconstruction  of  a  point  source  at  (0.25A.15A)  with  a  barrier  extending  from 
10A);  contour  lines  are  at  25%,  50%,  and  75%  of  the  maximum. 
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Figure  4.31.  Mismatch  reconstruction  of  a  point  source  at  (0.25A.15A)  with  a  barrier  extend 
(0, 10A)  to  (-°°,10A).  Values  below  50%  of  the  maximum  have  been  removed  for  clarity. 


lismalch  reconstruction  of  a  point  source  at  (— 2.25A.15A)  with  a  barrier  extending  from 
10A);  contour  lines  are  at  25%,  50%,  and  75%  of  the  maximum. 


Figure  4.33.  Mismatch  reconstruction  of  a  point  source  at  (-5.25A.15A)  with  a  barrier  extending  from 
(0,10A)  to  (-oo.lOA);  contour  lines  are  at  25%,  50%,  and  75%  of  the  maximum. 
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Figure  4.34.  Pattern-inalch  reconstruction,  using  a  false  barrier  at  Z=11A,  of  a  point  source  at 
(0.25A.15A)  with  a  true  barrier  at  Z=10A;  contour  lines  are  at  10%,  20%,  ....  90%  of  the  maximum. 
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location,  at  Z  =  1 1 A .  A  comparison  with  the  corresponding  results  with  no 
error.  Figures  4.26  and  4.30,  indicate  that  the  error  did  not  significantly  affect 
the  reconstructions  (note  that  only  3  contour  levels  were  used  in  Figure  4.30  while 
9  levels  were  used  in  Figure  4.35). 

4.8  Locating  Two  Sources 

When  more  than  one  source  is  present,  a  linear  superposition  of  the  signals 
from  the  sources  is  recorded  at  each  receiver.  However,  since  it  is  difficult  to 
classify  the  pattern-match  and  mismatch  methods  as  either  linear  or  non-linear, 
it  is  not  clear  whether  multiple  sources  could  be  reconstructed. 

To  explore  this  question,  a  simulation  was  run  with  two  equal-strength 
and  in-phase  sources  located  at  (0.25A,  15A)  and  (10.25A,  15A).  The  result  using 
the  pattern-match  method  is  shown  in  Figure  4.36.  The  field  near  the  left  source, 
at  (0.25A,  15A),  is  nearly  identical  to  the  corresponding  area,  in  Figure  4.26,  when 
only  the  left  source  was  used.  However,  as  can  be  clearly  seen  in  Figure  4.37,  the 
two  sources  are  not  located  with  the  same  resolution  due  to  the  influence  of  the 
barrier. 


The  result  using  the  mismatch  method  is  shown  in  Figure  4.3S.  Again,  the 
field  near  the  left  source  is  nearly  identical  to  the  corresponding  area,  in  Figure 
4.30,  when  only  the  left  source  was  used  (note  that  only  3  contour  levels  were 
used  in  Figure  4.30  while  9  levels  were  used  in  Figure  4.3S).  Figure  4.39  shows 
that  while  both  sources  have  been  located  with  nearly  the  same  resolution,  their 
apparent  strengths  are  not  equal. 


match  reconstruction  of  2  sources,  at  (0.25A.15A)  and  (10.25A.15A),  with  a  barrier 
to  (-°°,10AV,  contour  lines  are  at  10%,  20%,  ....  90%  of  the  maximum. 


Figure  4.37.  Pattern— match  reconstruction  of  2  sources,  at  (0.25A.I5A)  and  ( 10.25A.15A),  with  a  barrier 
extending  from  (0,1 0A )  to  (-°°,10A). 
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Figure  4.39.  Mismatch  reconstruction  of  2  sources,  at  (0.25A.15A)  and  (10.25A.15A),  with  a  barrier 
extending  from  (0.10A)  to  (— 10A) .  Values  below  50%  of  the  maximum  have  been  removed  for  clari 
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4.9  Execution  Times 

The  time  needed  to  complete  the  reconstruction  was  measured  using  the 
computer  system's  internal  clock,  which  is  claimed  to  be  accurate  to  within  0.01 
seconds.  Note  that  the  measurements  are  made  so  that  input,  output,  plotting, 
and  overall  normalization  operations  are  excluded.  However,  the  time  needed  to 
perform  the  normalization  of  the  fictitious  source  signals  for  the  pattern-match 
method  is  included.  The  reconstructions  were  performed  at  a  varying  number 
of  points,  depending  upon  the  situation;  therefore,  the  total  time  for  each  re¬ 
construction  has  been  divided  by  the  number  of  points  so  that  the  times  can 
be  directly  compared.  The  results  are  shown  in  Table  4.2  and  can  be  used  to 
evaluate  the  relative  speed  of  the  various  methods. 

The  speed  of  the  holography  method  is  the  same  when  the  source  is  par¬ 
tially  or  fully  hidden,  or  when  it  is  yields  a  reflection,  since  it  does  not  vary 
its  calculations  as  the  environment  varies.  The  pattern-matching  methods  are, 
as  expected,  significantly  slower  than  holography  since  they  cannot  take  advan¬ 
tage  of  the  efficiencies  of  the  FFT.  However,  there  is  no  obvious  reason  why  the 
mismatch  method  was  generally  faster  than  the  pattern-match  method  —  the 
mismatch  algorithm  is  slightly  more  complicated  and  so  it  was  expected  that  it 
would  be  slightly  slower. 

These  times  are  for  execution  on  the  MicroVAX  II  processor  and  floating 
point  accelerator  manufactured  by  the  Digital  Equipment  Corporation  (DEC). 
If  a  modern  super  computer  or  array  processor  had  been  used  instead,  the  time 
to  reconstruct  an  image  could  probably  be  reduced  by  at  least  two  orders  of 
magnitude. 
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Table  4.2.  Time  to  reconstruct  the  images. 


Environment 

Source 

Method 

Total  Time 

Time/Point 

free 

(0.25A,  15A) 

holography 

0.2m 

2.5ms 

free 

(0.25A.15A) 

pattern  match 

3.1m 

37.7ms 

free 

(0.25A,  15A) 

mismatch 

3.1m 

37.2ms 

reflections 

(— 5.25A,  7.5A) 

holography 

0.2m 

2.5ms 

reflections 

(— 5.25A,  7.5A) 

pattern  match 

24.1m 

-- 

396. Sms 

reflections 

(— 5.25A,  7.5A) 

mismatch 

24.1m 

395.9ms 

hidden 

(0.25A,  15A) 

holography 

0.2m 

2.5ms 

hidden 

(0.25A,  15A) 

pattern  match 

32.0m 

390.2ms 

hidden 

(0.25A,  15A) 

mismatch 

31.4m 

3S2.5ms 

hidden 

(-2.25A,  15A) 

holography 

0.2m 

2.5ms 

hidden 

(— 2.25A,  15A) 

pattern  match 

31.7m 

3S6.0ms 

hidden 

(— 2.25A,  15A) 

mismatch 

31.4m 

3S2.9ms 

hidden 

(— 5.25A,  15A) 

holography 

0.2m 

2.5ms 

hidden 

(— 5.25A,  15A) 

pattern  match 

31.4m 

382.2ms 

hidden 

(-5.25A,  15A) 

mismatch 

31.5m 

3S2.9ms 

A  -  wavelengths  m  -  minutes  ms  -  milliseconds 
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Chapter  5 

Swept-Frequency  Imaging 

One  of  the  purposes  of  this  study  was  the  development  of  a  new  method 
of  imaging  which  is  based  upon  a  sweep  of  the  transmitted  signal  through  a 
series  of  frequencies.  This  chapter  develops  the  method  and  conducts  a  variety 
of  simulations  in  order  to  explore  its  capabilities. 

5.1  Motivation  for  Sweeping  the  Frequency 

A  number  of  the  methods  examined  by  this  study  make  use  of  multiple 
frequencies.  When  there  are  several  sources  differing  in  their  spectral  content,  the 
eigenvector  approach  uses  this  fact  to  help  separate  the  information  belonging 
to  the  sources  and  thereby  increase  its  resolving  capabilities.  The  correlation 
method  uses  pairs  of  receivers  to  estimate  the  relative  time  delay  and  bearing, 
and  so  more  precisely  locates  the  peak  in  the  correlation  as  the  bandwidth  is 
increased.  Other  methods,  such  as  beamforming  and  holography,  can  process 
only  one  frequency  at  a  time,  and  so  can  only  combine  the  results  for  multiple 
frequencies  at  the  end.23,24 

One  thing  which  all  these  approaches  have  in  common  is  that  they  con¬ 
centrate  on  the  signal  processing  viewpoint,  treating  the  information  at  each 
frequency  as  being  unrelated  to  the  information  at  other  frequencies.  However, 
if  we  look  at  the  underlying  physics,  it  is  obvious  that  such  a  situation  would  be 
unusual,  if  not  impossible,  for  scatterers  of  acoustic  waves  since  their  scattering 
must  depend  upon  the  frequency  in  some  deterministic  if  unknown  manner. 


Such  a  frequency  dependence  was  exploited  by  Chan  et  al to  emulate 
motion  by  the  receivers  along  “equivalent  scan  vectors”  so  that  the  effective  size 
of  the  receiving  array  was  increased.  A  modest  body  of  work  was  produced  by  this 
group,'6,77’'8  along  with  observations  by  other  parties  that  some  of  the  reported 
results  were  unacceptable'9  or  exhibited  ambiguities.80  Although  a  lack  of  further 
publications  seems  to  indicate  that  this  work  has  not  been  pursued,  it  seems  clear 
that  in  some  situations  the  effective  size  of  the  array  in  one  of  the  dimensions  can 
be  increased  by  a  factor  proportional  to  the  bandwidth  of  the  frequency  sweep. 

One  of  the  weaknesses  of  Chan’s  approach  is  that  the  calculation  of  the 
equivalent  scan  vectors  requires  that  the  center  of  the  scattering  object  be  known. 
Figures  5.1  and  5.2  depict  the  coordinate  system  and  scan  vectors  for  a  rectan¬ 
gular  scatterer  with  its  center  on  the  Z  axis.  It  is  obvious  that  the  scan  vectors 
cannot  be  calculated  if  the  center  of  the  scatterer  is  not  known. 

A  solution  to  this  limitation  presents  itself  when  we  consider  the  abilities  of 
the  pattern-match  method  discussed  in  Chapters  2  and  3,  in  place  of  the  concept 
of  scan  vectors. 

In  the  pattern-matching  method,  rather  than  assume  that  we  know  the 
center  of  the  scatterer,  we  move  the  fictitious  scatterer  through  all  locations  in 
object  space,  i.e.,  all  locations  which  could  possibly  contain  an  actual  scatterer. 
At  each  of  these  locations,  the  received  signals  are  compared  to  what  would  have 
been  received  from  the  fictitious  scatterer,  at  the  given  frequency. 

However,  to  the  extent  that  the  frequency-dependence  of  the  scatterer  is 
predictable,  another  type  of  match  could  be  performed.  Using  the  received  pat- 
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Definition  of  the  equivalent  scan  vectors  in  a  proposed  swept-frequency  imaging  system. 
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tern  across  the  array  at,  say,  the  lowest  frequency,  we  can  predict  how  the  pattern 
will  shift  and  scale  as  the  frequency  is  swept  up  to  higher  frequencies,  given  the 
location  of  the  fictitious  scatterer  and  the  known  location  of  the  transmitter  and 
receivers.  The  degree  to  which  the  actual  and  predicted  signals  match,  summed 
over  the  entire  frequency  sweep,  then  constitute  the  total  match  at  that  location 
in  object  space.  As  before,  we  then  interpret  the  best  match  or  matches  in  object 
space  as  the  location  of  the  source  or  sources. 

5.2  Scattering  Models 

The  purpose  of  the  scattering  model  used  in  this  chapter  of  the  study  is  to 
provide  the  frequency-swept  imaging  simulations  with  “received”  signals  which 
should  mimic  features  of  real,  physical  phenomena  with  reasonable  accuracy  when 
these  features  may  have  an  effect  on  the  imaging  method.  Thus,  the  frequency 
dependence  of  the  scattering  should  be  faithfully  simulated.  On  the  other  hand, 
the  absolute  levels  of  the  simulated  signals  is  probably  irrelevant. 

It  seems  that  two  categories  of  frequency-dependent  scattering  would  suf¬ 
fice.  First,  a  pair  of  closely  spaced  point  scatterers  will  be  used  to  represent 
specular  reflections.  While  an  individual  point  scatterer  has  no  frequency  de¬ 
pendence,  a  system  of  two  or  more  will  exhibit  frequency-dependent  interference 
patterns.  Second,  a  rigid  plate  reflector  will  be  used  to  represent  an  extended 
scatterer. 

5.2.1  Specular  Scatterers 

The  expression  for  point  scatterers  is  very  nearly  identical  to  the  general 
case  given  as  part  of  the  definition  of  terms  in  Eq.  (1.2).  For  the  following  sim- 
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ulations,  we  simplify  the  formula  by  setting  the  scattering  strength  to  unity  so 
that  Sm  ~  1-  The  transmitted  signal  is  chosen  to  be  a  continuous  sine  wave  of 
unity  amplitude  so  that  the  time-shifted  signal,  E  (ti  —  T^j ,  can  be  replaced 
by  a  phase-shifted  signal,  expf — j  '2n  f  —  T-^Jj  with  Tfm  given  by  Eq.  (1.5). 
Finally,  since  the  time  at  which  we  sample  the  complex  received  signal  is  irrele¬ 
vant  while  using  continuous  sine  waves  (as  long  as  all  the  receivers  are  sampled 
at  the  same  time),  we  can  set  t/  =  0  so  that  the  expression  for  the  received  signal 
becomes 


R,  (A'f )  =  £ 


^  exp  (  — ?  k  ( Dim  +  DTm )) 


m=l 


Dxm  ’  Dj'f, 


(5.1) 


The  distance  Dim  is  that  between  the  mth  reflector  and  the  1th  receiver  located 
at  (X-1,  Z?)  where  we  have  chosen  to  locate  our  linear  receiver  array  along  the 
Z-axis  so  that  Zf1  =  0.  The  distance  D^m  is  that  between  the  mth  reflector 
and  the  point  source  at  (XT ,ZT^j.  These  coordinates  and  distances  are  shown 
in  Figure  1.1,  and  the  distances  are  given  by  Eq.  (1.3)  and  Eq.  (1.4). 


Eq.  (5.1)  was  implemented  in  the  FORTRAN  subroutine  and  then  used  to 
generate  Figures  5.3  and  5.4  for  the  case  of  two  reflectors  separated  by  12".  In 
these  figures,  the  reflectors  were  centered  around  the  origin  along  the  X-axis, 
and  the  reflected  field  was  evaluated  at  a  radius  of  100"  from  the  origin.  Several 
other  source  locations  were  also  examined  but  not  shown  here;  in  all  of  these 
other  cases,  the  magnitude  and  phase  peaks  moved  as  expected  from  0°  to  the 
negative  of  the  source  angle. 

5.2.2  Extended  Scatterers 


The  classical  theories  of  Huygens  shall  be  used  to  construct  a  simulation 
of  the  scattering  of  an  acoustic  wave  by  a  finite  rigid  plate.  The  incident  wave  is 


source 


12S 


the  origin  and  separated  by  12".  The  wavelength  is 


modeled  as  a  harmonic,  spherically  spreading  one  originating  at  a  point  source. 
We  are  interested  in  the  field  reflected  back  towards  our  linear  receiver  array 
which  is  in  the  general  vicinity  of  the  source;  that  is,  the  model  does  not  need 
to  consider  diffraction  around  the  plate.  Thus,  we  can  apply  Huygen’s  principle 
wherein  each  point  on  the  plate  is  considered  to  be  a  point-source  of  radiation 
with  a  phase  and  amplitude  equal  to  that  of  the  incident  wave. 

To  be  rigorous,  the  signal  at  each  receiver  should  not  only  include  waves 
reflected  from  the  plate  but  also  waves  arriving  directly  from  the  transmitter. 
We  shall  delete  this  direct  path  in  this  simulation.  In  an  experiment,  the  signals 
arriving  along  the  direct  path  could  be  subtracted  out  since  the  exact  location  of 
the  transmitter  and  receivers,  and  the  speed  of  sound,  is  known.  Or,  windowed 
sine  waves  could  be  used  instead  of  a  continuous  wave  so  that  the  two  signals  can 
be  distinguished  by  a  time  window  at  the  receivers. 

After  a  presentation  of  the  integral  and  discrete  expressions  for  the  scat¬ 
tered  waves,  subsequerd  sections  will  exercise  the  simulation.  An  extra  amount 
of  effort  was  put  into  testing  the  algorithm  as  part  of  an  attempt  to  explain  an 
unexpected  scattering  pattern  which  occurs  when  the  plate’s  width  becomes  large 
relative  to  a  wavelength. 

5. 2. 2.1  Scattering  from  a  Rigid  Plate 


The  Huygens-Rayleigh  integral  for  scattering  from  a  plate  can  be  expressed 
in  terms  of  the  incident  pressure81  as 

Xmax  Ymax 


R 


I  I 


X  Y 
^ min  1  min 


Pl(x,y)  exp  {-j  k  £>,•) 

Dx 


dydx.  (5.2) 


where  the  integrals  are  over  the  surface  of  the  plate.  In  this  study,  we  are  con¬ 
sidering  only  the  two  dimensional  X-Z  plane  and  so  the  dependence  and  integral 
over  Y  will  be  dropped.  Also,  we  have  chosen  to  locate  our  linear  receiver  array 
along  the  Z-axis  and  so  Z '/*  =  0.  The  incident  pressure  Pj  is  given  by 


Pl(x)  =  P0 


exp  (-j  k  Dt ) 


The  distances  Dt  and  Dt  are  given  by  Eq.  (1.3)  and  Eq.  (1.4),  respectively,  except 
that  the  location  of  the  sources,  is  now  x,  the  variable  of  integration  —  i.e., 
the  location  of  Huygen’s  equivalent  sources  on  the  plate.  With  these  assumptions 
and  definitions,  Eq.  (5.2)  becomes 


There  does  not  seem  to  be  an  exact  closed-form  solution  to  this  integral. 
Approximations  to  the  square  root  in  the  D{  and  Dt  can  allow  the-jntegral  to 
be  evaluated  in  the  Fraunhofer  zone,  or  when  the  receiver  is  near  the  axis  of  the 
plate82.  Unfortunately,  the  resulting  expressions  would  be  invalid  for  some  of 
the  arrangements  we  wish  to  consider  for  the  frequency-swept  imaging  simula¬ 
tions.  Thus,  we  shall  evaluate  the  integral  numerically  with  a  simple  but  robust 
summation.  If  we  wish  to  break  the  integral  into  NP  divisions,  we  let 

dx  —y  AX, 


AX  = 


A  max  Amm 
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Then,  Eq.  (5.4)  becomes 

n/,'R\  jkPat\X  exp(-jk(Dim  +  Dm,)) 

R‘  E- )  ~  —oT-  s.,  •  cTm 

where 

Xi  =  xmtn  +  AX  •  (m  -  0.5)  (5.7) 

This  summation  is  arranged  so  the  the  integral  over  each  division  is  approximated 
by  the  value  of  the  integral  at  the  center  of  the  division. 


5. 2. 2. 2  Plate  Scattering  Results 


The  simulation  of  plate  scattering  was  first  exercised  for  the  case  of  a  12” 
plate  located  at  the  origin.  In  Figures  5.5  and  5.6,  the  magnitude  and  phase, 
respectively,  of  the  reflected  pressure  are  shown  due  to  a  source  100"  from  the 
origin  at  0°.  Figure  5.7  shows  a  similar  case  except  that  the  source  is  at  —45°.  In 
all  of  these,  the  field  was  evaluated  at  locations  which  were  100"  from  the  origin, 
at  angular  increments  of  0.25°  —  that  is,  each  plot  consists  of  721  points. 

5. 2. 2. 3  Comparison  to  Closed-Form  Expressions 

It  would  be  encouraging  to  be  able  to  compare  the  figures  of  the  previous 
section  to  the  closed-form  expressions  obtained  by  making  the  Fraunhofer  approx¬ 
imation.  This  approximation  is  made  in  the  distances  D{  and  Dt  in  Eq.  (5.4). 
A  small  error  in  the  distances  will  cause  only  a  small  error  where  they  are  used 
in  the  denominator  of  Eq.  (5.4).  On  the  other  hand,  the  complex  exponential  is 
very  sensitive  to  such  errors;  and,  the  problem  is  made  worse  by  the  multiplicative 
factor  k. 
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summation  used  200  terms  of  width  0.06"  (0.06A). 


source 


Figure  5.6.  Phase  of  the  pressure  from  a  13,045.2  Hz  source  at  0°,  scattered  by  a  12  plate.  The 
summation  used  200  terms  of  width  0.06"  (0.06X). 
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Figure  5.7.  Magnitude  of  the  pressure  from  a  13,045.2  Hz  source  at  -45°,  scattered  by  a  12  plate.  The 
summation  used  200  terms  of  width  0.06"  (0.06A). 
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If  we  wish  to  limit  the  error  in  the  phase  of  this  exponential  to,  say,  10°, 
we  must  limit  the  error  in  the  distances  to 


,  r  radians  1  „ _ 

10°  •  — — - T  =  0.02S 

180°  k 


(5.8) 


where  we  have  set  A  =  1"  so  that  k  =  2tt.  We  are  interested  in  simulations 
of  a  12"  plate;  thus,  the  greatest  error  will  occur  when  computing  the  distance 
between  the  edge  of  the  plate  and  a  field  point  at  some  distance  Z  at  some  angle  0. 
This  distance  can  be  computed  without  and  with  the  Fraunhofer  approximation 


as 


D exact  -  \J  (z  COS (0)  -  b)  “  +  Z2 , 

D Fraunhofer  —  Z  6  sin(0). 


(5.9) 

(5.10) 


These  two  distances  were  compared  for  various  values  of  Z  as  0  was  increased 
to  discover  the  maximum  angle  for  which  the  Fraunhofer  approximation  met  the 
error  limit.  The  results,  shown  in  Table  5.1,  indicate  that  around  and  below 
Z  =  600",  the  Fraunhofer  approximation  fails  the  criteria  at  all  angles.  Since 
we  are  interested  in  using  this  simulation  for  distances  of  around  Z  =  100", 
comparison  to  closed-form  expressions  using  the  approximation  cannot  be  used. 

5. 2. 2. 4  Numerical  Suitability  of  the  Summation 

One  of  the  common  causes  of  a  loss  of  numerical  accuracy  is  that  large 
and  small  numbers  are  added  together  in  an  inappropriate  order,  so  that  the 
contributions  from  the  small  numbers  are  lost.  In  order  to  explore  this  possibility, 
the  scattering  shown  in  Figure  5.5  was  recalculated  with  the  accumulated  complex 
value  of  the  summation  shown  in  Figure  5.8,  and  the  individual  complex  terms 
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Table  5.1.  Maximum  angle  at  which  the  Fraunhofer  approximation  is  valid. 


z 

o 

<600.0  inches 

none 

700.0  inches 

0.1350  degrees 

1,000.0  inches 

0.2521  degrees 

10,000.0  inches 

0.1310  degrees 

100,000.0  inches 

0.037S  degrees 

1,000,000.0  inches 

0.0184  degrees 

of  the  summation  plotted  in  Figures  5.9.  Note  that  the  width  of  the  plate  was 
increased,  from  that  used  in  earlier  figures  to  24",  so  as  to  better  illustrate  the 
process. 

The  summation  shown  in  Figure  5.8  begins  at  zero;  therefore,  the  first 
point  on  the  curve  (at  the  end  with  the  heavy  line)  is  the  same  as  the  first  term. 
The  last  point  on  this  curve  is  the  final  value  of  the  summation.  In  Figure  5.9, 
only  half  the  terms  seem  to  be  present.  Actually,  they  are  all  present:  because  of 
symmetry,  the  terms  for  one  half  of  the  plate  are  exactly  equal  to  the  terms  for 
the  other  half,  and  so  pairs  of  points  plot  exactly  on  top  of  one  another. 

These  figures  indicate  that  the  evaluation  of  Eq.  (5.6)  should  be  a  numer¬ 
ically  accurate  process  since  the  summation  terms  have  nearly  the  same  magni¬ 
tude,  and  are  within  approximately  one  order  of  magnitude  of  the  final  result. 


137 


j,OI  ♦  lejjfeiui  jo  jjbj  /(jbuiSbuii 


Figure  5.8.  Accumulated  Huygens-Rayleigh  sum  as  we  move  across  the  plate  from  X=-12"  (the  end  of 
the  curve  with  the  heavy  line)  to  X=+12”.  The  receiver  and  13,045.2  Hz  source  were  100”  away  at  0°. 
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5. 2. 2. 5  Number  of  Terms 

The  accuracy  by  which  the  summation  of  Eq.  (5.6)  approximates  the  in¬ 
tegral  of  Eq.  (5.4)  is  governed  by  the  number  of  terms  Np-,  therefore,  it  seems 
prudent  to  examine  the  convergence  of  the  sum  as  Np  is  increased.  As  illustrated 
in  Figures  5.S  and  5.9,  the  summation  is  essentially  one  of  vector  addition  in  the 
complex  plane  where  the  length  of  the  vectors  are  nearly  equal.  However,  the 
direction  of  these  vectors  will  vary  more  rapidly  as  the  frequency  or  distances  are 
increased. 

A  series  of  tests  were  run  where  a  12 "-wide  plate  scattered  waves  from  a 
13,045.2  Hz  source  located  on  the  axis  of  the  plate  at  a  distance  of  100".  The 
normalized  field  at  a  point  100"  away  at  an  angle  of  45°  is  reported  in  Table 
5.2.  The  results  obtained  using  10,000  terms  was  used  as  the  reference  value  for 
normalization,  which  was  done  separately  for  the  real  and  imaginary  parts.  Given 
that  the  wavelength  is  1",  the  results  for  using  only  a  small  number  of  divisions 
are  remarkably  good. 


Table  5.2.  Normalized  complex  signal  at  13,045.2  Hz  and  45°. 


Np 

Width  of  Term 

Real  Part 

Imaginary  Part 

30  terms 

0.4000  inches 

1.1410 

1.1271 

60  terms 

0.2000  inches 

1.0327 

1.0301 

100  terms 

0.1200  inches 

1.0116 

1.0100 

1,000  terms 

0.0120  inches 

1.0002 

1.0000 

10,000  terms 

0.0012  inches 

1.0000 

1.0000 

HO 


A  similar  series  of  test  were  then  run  where  the  field  point  was  moved  to 
0°  and  the  frequency  was  increased  to  30,000  Hz  so  that  the  wavelength  is  now 
0.435”.  The  results,  shown  in  Table  5.3,  also  show  good  convergence.  Because  of 
these  tests,  it  was  decided  to  set  Np  =  200  for  the  simulation  of  plate  reflection 
in  this  study  although  higher  numbers  might  be  used  at  very  high  frequencies. 

Table  5.3.  Normalized  complex  signal  at  30,000  Hz  and  0°. 


Np 

Width  of  Term 

Real  Part 

Imaginary  Part 

10  terms 

1.2000  inches 

0.9926 

0.9255 

50  terms 

0.2400  inches 

0.9995 

0.9974 

100  terms 

0.1200  inches 

0.9999 

0.9993 

1,000  terms 

0.0120  inches 

1.0000 

1.0000 

10,000  terms 

0.0012  inches 

1.0000 

1.0000 

5. 2. 2. 6  Results  Using  Higher  Precision 

The  simulation  of  reflections  from  a  plate  was  implemented  in  “double¬ 
precision”  arithmetic  which  means  that  the  computer’s  internal  memory  repre¬ 
sents  each  number,  or  each  part  of  a  complex  number,  with  an  accuracy  equal 
to  approximately  16  decimal  digits.  Maintenance  of  this  accuracy,  of  course, 
depends  upon  the  numerical  methods  employed  in  an  algorithm.  However,  the 
support  routines,  such  as  for  the  evaluation  of  the  sine  and  cosine  functions,  are 
designed  to  deliver  this  accuracy. 


The  computer  being  used  can  also  support  “quadruple-precision”  arith- 
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metic  where  each  number  and  the  support  routines  yield  approximately  33  deci¬ 
mal  digits  of  accuracy.  A  special  version  of  the  simulation  was  prepared  and  the 
values  shown  in  Table  5.3  were  recalculated  using  this  higher  precision.  A  com¬ 
parison  between  the  normalized  quadruple-  and  double-precision  results  revealed 
no  differences  within  the  first  5  digits.  Because  of  these  tests,  it  was  decided  that 
the  faster  double-precision  arithmetic  was  adequate. 

5. 2. 2. 7  Scattering  at  Large  Width- Wavelength  Ratios 

In  preparation  for  the  imaging  simulations,  the  plate  simulation  was  ex¬ 
ercised  by  varying  the  plate  width,  location,  and  the  frequency  of  the  source. 
At  higher  frequencies,  an  anomaly  occurred.  We  expected  the  scattering  pattern 
shown  in  Figure  5.5  to  be  essentially  reproduced  at  higher  frequencies  except 
that  the  pattern  would  be  compressed  towards  the  angle  of  incidence,  which  was 
0°  in  these  exercises.  However,  Figure  5.10  shows  that  we  instead  obtained  a 
pattern  with  a  dip  where  we  expected  the  central  peak.  To  explore  this  further, 
the  pattern  across  an  array  of  receivers  was  generated  for  a  sweep  of  frequencies. 
The  results,  shown  in  Figure  5.11,  indicate  that  while  the  sidelobes  do  indeed 
compress  as  expected  as  the  frequency  is  increased,  the  anomalous  region  not 
only  continues  but  increases  in  size. 

Other  tests  had  revealed  that  the  width  of  the  plate  was  also  a  factor  in 
these  unexpected  results.  In  Figure  5.12,  we  have  reproduced  the  simulation  used 
in  Figure  5.5  except  that  the  width  of  the  plate  has  been  increased  from  12"  to 
48".  In  order  to  observe  the  progressive  effect  of  plate  width,  the  width  was  swept 
from  12"  to  4S",  in  steps  of  0.5",  and  the  normalized  results  displayed  in  Figure 
5.13.  Note  that  as  the  width  of  the  plate  is  varied,  the  width  of  the  anomalous 


source 
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ation  used  200  terms  of  width  0.06"  (0.14X). 


re  5.11.  Magnitude  of  the  pressure  from  a  source  at  the  origin,  scattered  by  a  12"  plate  100"  away, 
summation  used  200  terms  of  width  0.06".  The  wavelength  varies  from  1.305"  to  0.177". 


source 


lion  used  200  terms  of  width  0.24"  (0.24A). 


summation 
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region  varies  proportionally. 

We  can  gain  a  better  understanding  of  these  simulations  by  applying  Babi- 
net’s  principle;83  that  is,  we  replace  the  rigid  plate  in  its  plane  with  an  infinite 
baffle  in  the  same  plane  which  contains  an  aperture  where  the  plate  had  formerly 
been  located.  In  addition,  we  replace  the  point-source  transmitter  with  its  mirror 
image,  located  on  the  opposite  side  of  the  plane.  In  this  alternative  arrangement, 
the  phase  and  amplitude  at  each  location  in  the  new  aperture  is  identical  with 
that  which  had  formerly  existed  at  the  same  locations  on  the  surface  of  the  plate. 
Thus,  Eq.  (5.4)  continues  to  describe  the  signal  at  the  receivers. 

The  effects  seen  in  Figure  5.13  can  now  be  interpreted  as  the  field  of  a  point 
source  after  passing  through,  and  experiencing  diffraction  at,  an  aperture.  For 
the  special  case  of  the  point  source  location  being  on  the  axis  of  the  aperture  (i.e., 
opposite  its  center),  Skudrzyk84  noted  that  the  field  consists  of  the  interference 
between  the  incident  wave  and  a  contribution  that  originates  at  the  boundary 
of  the  aperture.  The  special  case  applies  here  since  the  transmitter  has  been 
located  on  the  Z-axis  about  which  the  plate  has  also  been  centered.  Given  the 
speed  of  sound  c  =  13045.2  inches/second,  we  find  that  the  edge  of  the  first  or 
central  Huygens  zone  to  be  10.01,  6.543,  or  4.214  inches  from  the  origin  when 
the  frequency  is  13045.2,  30,500,  or  73,500  cycles/second,  in  that  order.  This 
corresponds  to  zones  with  a  diameter  (or  width,  in  our  two-dimensional  case)  of 
20.02,  13.0S6,  or  8.428  inches,  in  that  order. 

It  therefore  seems  reasonable  that  the  earlier  examples,  such  as  Figure 
5.5,  did  not  exhibit  the  interference  between  the  incident  and  boundary  waves 
because  the  plate  was  not  large  enough  to  encompass  even  the  first  Huygens  zone. 
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The  small  dip  where  the  central  peak  had  been  expected  in  Figure  5.10  is  not  an 
anomaly  but  is  instead  due  to  the  plate  between  just  wide  enough,  relative  to  a 
wavelength,  for  this  interference  to  be  observed. 

5.3  Interpolation 

The  central  idea  in  the  swept-frequencv  method  developed  in  the  next 
section  is  the  comparison  of  the  pattern  recorded  across  the  receiving  array  at 
one  frequency  to  that  recorded  at  another  frequency.  The  approach  (discussed 
later)  shall  be  to  use  the  higher  of  the  two  frequencies  directly.  Thus,  we  shall 
need  to  compress  the  pattern  at  the  lower  frequency  to  estimate  how  it  would  look 
at  the  higher  frequency.  That  is,  the  signal  Ri  recorded  at  the  lower  frequency 
at  location  X-1  in  the  linear  array  is  moved  to  a  new  location  Xj1  representing 
the  compression  caused  by  the  increase  in  frequency. 

In  general,  the  new  locations  „Y,-  will  not  happen  to  be  equal  to  the  original 
receiver  locations  at  which  the  higher  frequency  is  known.  We  will  therefore  need 
to  interpolate  on  the  compressed  pattern  to  find  its  value  at  the  original  receiver 
locations.  While  this  is  in  some  ways  a  minor  point  in  the  development,  it  is 
central  to  the  analysis  not  only  in  terms  of  the  accuracy  of  the  results  but  also 
in  terms  of  the  speed  of  execution. 

A  piecewise-cubic  Hermite  interpolation  algorithm  was  adapted  from  an 
algorithm  given  by  Conte  and  de  Boor85  and  implemented  as  a  fortran  subroutine. 
This  algorithm  requires  the  slope  of  the  function  at  :h  point.  Since  the  slope 
of  the  scattering  pattern  across  the  array  cannot  be  ctly  measured,  the  slope 
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y  •  is  estimated  as 

/  yi+ 1  -  yi- 1 

y*  =  - - 

•£»+i  3-t—i 

where  the  yi  are  the  function’s  values  at  locations  a:,-. 

It  was  felt  that  it  would  be  better  for  this  operation  to  return  no  value 
than  to  blindly  estimate  answers  which  were  not  robustly  supported  by  the  given 
data;  therefore,  the  interpolation  subroutine  signals  that  it  has  set  the  result  to 
zero  if  the  requested  point  of  interpolation  is  outside  certain  limits.  The  limits 
have  been  chosen  so  that  there  is  always  at  least  two  points  between  the  point  of 
interpolation  and  the  beginning  or  end  of  the  table  of  given  data.  These  limits 
were  selected  because  the  algorithm  needs  one  point  on  either  side  of  the  point 
of  interpolation  for  the  calculation  of  slope. 

Since  the  pattern  across  the  array  often  looks  like  a  sine  wave,  a  test  was 
run  using  9  samples  of  a  sine  wave  every  45°.  Interpolated  values,  shown  in 
Figure  5.14,  were  requested  every  1°.  As  expected,  zeroes  were  returned  beyond 
the  interpolation  limits.  The  maximum  error,  occurring  near  the  peaks,  are 
acceptable  for  this  study. 

5.4  Development  of  the  Analysis  Method 

In  Section  5.2,  the  scattering  patterns  from  point  and  extended  scatterer s 
were  examined  for  correctness;  now,  we  examine  them  again  to  determine  their 
frequency  dependence  which  will  then  be  used  to  develop  the  new  swept-frequency 
imaging  method. 
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Figure  5.14.  Values  interpolated  using  a  piecewise  cubic  polynomial  (line)  and  the  control  points  (large 
circles)  used  in  the  interpolation.  Additional  values  (small  circles)  are  shown  for  comparison. 
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5.4.1  Frequency  Dependence 

We  wish  to  develop  a  model  of  the  frequency  dependence  of  a  scattering 
pattern.  The  approach  shall  be  to  develop  the  model  for  the  rather  simple  case 
of  two  point  reflectors,  but  then  test  the  model  for  both  pairs  of  point  reflectors 
as  well  as  plates. 


The  basic  geometry  for  determining  the  signal  scattered  by  a  point  reflec¬ 
tor  is  given  in  Figure  1.1.  The  signal  at  a  particular  receiver  located  at  (XR,  ZR) 
due  to  one  reflector  is  given  by  Eqs.(1.2)  to  (1.5)  with  Ns  =  2.  For  a  single¬ 
frequency  source,  the  time  delays  can  be  written  as  phase  delays.  Then,  we 
merely  sum  the  contributions  from  the  two  reflectors  so  that  we  have 


R{Xr,  Zr ) 


exp (— j  k  (Pti  +  Dr\))  exp(-j  k  (DT2  +  Dr 2)) 
Dt\  •  Dr  1  Dt2  ■  DR2 


(5.12) 


where  the  Drs  are  the  distances  between  the  transmitter  and  scatterer  s  and  the 
Drs  are  the  distances  between  scatterer  s  and  the  receiver.  When  the  distance 
d  between  the  scatterers  is  small  relative  to  the  distance  to  the  source  and  to 
the  receiver,  we  can  make  an  approximation  where  we  treat  the  pressure  waves 
as  being  planar  instead  of  spherical.  This  lets  us  write  the  previous  distances  in 
terms  of  the  distances  to  the  point  midway  between  the  reflectors.  If  Dr  is  the 
distance  to  this  point  from  the  transmitter,  and  Dr  is  the  distance  to  this  point 
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from  the  receiver,  we  can  write 

Dti  — »  Dt  —  T-dsin  i 

Dj2  —>  Dt  +  \d  sin  (d>T)  , 

(5.13) 

Dr\  Dr-  ^dsin  , 

Dr2  — >  Dr  +  \dsm  (<f>R'j 

where  4>T  and  <frR  are  the  angles  to  the  transmitter  and  receiver,  respectively, 
relative  to  the  normal  to  the  line  connecting  the  reflectors.  Implicit  in  these 
equations  is  the  fact  that  the  line  connecting  the  point  reflectors  is  parallel  to  the 
X-axis,  and  that  reflector  2  lies  at  an  X-coordinate  greater  than  that  of  reflector 
1.  We  are  mainly  interested  in  the  scattering  pattern’s  angular  rather  than  ampli¬ 
tude  dependence  upon  frequency;  futhermore,  the  amplitude  terms  change  very 
slowing  with  frequency.  Therefore,  we  can  drop  the  amplitude  factors,  leaving 

R(Xr,Zr) 

~  exp  k  Dt  +  Dr  -  \d  ^sin  ( <f>T )  +  sin 

+  exp  j  k  ^ Dt  +  Dr  +  kd  ^sin  (cfi1"}  +  sin  (cj>R , 

=  2  exp(— j  k  (Dt  +  Dr)^  •  cos  ^ k  d  ^sin  +  sin 

The  term  cos(. . .)  in  Eq.  (5.14)  describes  the  directivity  of  the  scattering  pattern 
at  a  given  frequency  /  where  k  —  We  wish  to  determine  how  the  directivity 
will  change  when  we  shift  to  a  new  frequency  f.  That  is,  if  some  particular 
feature  in  the  pattern  (a  peak,  a  null,  or  any  feature  inbetween)  is  at  angle  <j>R 
at  frequency  /,  to  what  new  angle  <j)R  does  the  feature  move  at  new  frequency 
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/?  It  is  convenient  to  answer  this  question  by  examining  the  nth  maximum  since 
the  entire  pattern  must  exhibit  the  same  frequency  dependence.  The  maximum 
points  in  the  pattern  will  occur  at 

U-d  (sin  (<^T)  +  sin  (<f>R)^j  =  n  tt  for  n  =  1,2,3,...,  (5.15) 


where  we  have  excluded  the  central  maximum  at  n  =  0.  The  parameters  n  and 
d  do  not  change  with  frequency.  Therefore,  we  solve  for  these  parameters  for 
frequencies  /  and  /  and  equate  the  results  so  that  we  have 

—  •  (sin  ( (f>T )  +  sin  =  ~  *  (sin  ( <t>T )  +  sin  (5.16) 


or 


=  arcsin^—  sin  (y4>T^j  +  y  '  (sin  ( <t>T )  + 


sin 


(5.17) 


If  the  scattering  center  (the  point  midway  between  the  two  reflectors)  is,  or  is 
assumed  to  be,  located  at  (XF ,  ZF),  then  we  can  convert  the  angular  frequency 
dependence,  given  by  Eq.  (5.17),  into  a  linear  frequency  dependence,  for  locations 
(XF,  O)  in  the  linear  receiver  array.  This  relationship  is 


XtR  =  XF  +  ZF  ■  tan  (^f) 


(5.18) 


where  <f>R  and  <f>R  are  the  angular  locations  of  the  z’th  receiver  relative  to  the 
scattering  center. 


The  ability  of  Eq.  (5.18)  to  represent  the  frequency  dependence  of  a  scat¬ 
tering  pattern  was  explored  by  first  calculating  the  pattern,  due  to  two  point 
reflectors,  at  a  linear,  equally  spaced  receiver  array  for  a  set  of  frequencies.  The 
results  are  shown  in  Figure  5.15  with  the  lower  90%  of  the  data  removed  so  that 


Patterns  at  an  array  on  the  X-axis  (64  receivers)  versus  frequency  (50  frequencies,  250  Hz 
3  2  point  reflectors  centered  at  (-20”, 100”)  and  a  source  at  (20”, 0). 
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the  location  of  the  peaks  can  be  more  easily  seen.  Then,  Eq.  (5. IS)  was  applied 
to  shift  all  frequencies,  except  30,000  Hz,  to  /  =  30,000  Hz.  The  shifted  results 
in  Figure  5.16  exhibit  two  phenomena.  First,  the  patterns  at  all  frequencies  are 
now  aligned,  indicating  that  Eq.  (5.18)  has  done  a  good  job  of  predicting  the 
frequency  dependence.  Second,  there  are  regions  of  the  shifted  pattern  which 
cannot  be  determined  (i.e.,  the  interpolation  returned  a  zero);  these  correspond 
to  regions  which  had  been  off  the  end  of  the  array  before  the  shift.  The  undeter¬ 
mined  region  is  greatest  where  the  shifting  is  greatest,  at  the  lower  frequencies. 
This  phenomenon  means  that  even  if  we  can  predict  the  frequency-induced  shift 
of  a  pattern,  the  results  may  be  less  than  useful  if  the  value  of  A /  (i.e.,  the 
difference  between  /  and  /)  is  too  great. 

The  same  information  shown  in  Figure  5.16  is  shown  again  in  Figure  5.17 
as  a  series  of  superimposed  curves.  A  repeating  cycle  of  solid,  long-dash,  and 
short-dash  lines  were  used  to  plot  the  curves,  beginning  at  the  highest  frequency. 
Note  that  only  the  highest  frequency,  a  solid  line,  extends  all  the  way  to  the  right- 
hand  side  of  the  plot  —  as  the  frequency  decreases,  the  undetermined  region  of 
the  shifted  pattern  increases  and  so  the  curve  ends  (goes  to  0)  at  some  decreasing 
location  in  X.  This  figure  exhibits  an  excellent  frequency  compensation  for  the 
lateral  location  of  the  pattern.  Unfortunately,  the  lack  of  any  frequency  compen¬ 
sation  for  the  amplitude  of  the  pattern  in  Eq.  (5.18)  can  be  seen  in  the  different 
levels  at  the  peaks  and  nulls. 

The  information  shown  in  Figure  5.17  has  been  regenerated  for  Figure 
5.18  except  that  the  pair  of  point  reflectors  has  now  been  replaced  with  a  plate 
12"  wide.  The  frequency  compensation  for  the  lateral  location  of  the  pattern 
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Figure  5.16.  Patterns,  due  to  2  point  reflectors,  at  an  array,  versus  frequency.  The  pattern  at  each 
frequency,  except  30,000  Hz,  has  been  shifting  so  as  to  approximate  a  pattern  at  30,000  Hz. 


15S 


is  remarkably  good,  especially  in  light  of  the  fact  that  Eq.  (5.18)  was  derived 
for  two  point  reflectors.  This  success  is  probably  due  to  the  fact  that  the  two 
edges  of  a  reflective  plate  are  the  major  contributors  to  its  scattering  pattern. 
The  frequency  dependence  of  the  amplitude  of  the  plate’s  pattern,  however,  has 
resulted  in  relatively  large  differences  at  the  peaks  in  the  pattern.  The  only  way 
in  which  this  difference  could  be  minimized  is  to  minimize  A/. 

5.4.2  Forming  a  Swept-Frequency  Image 

We  wish  that  the  input  to  this  new  imaging  method  be  a  set  of  data 
like  that  shown  in  Figure  5.15:  a  received  signal,  versus  receiver  location,  versus 
frequency.  We  also  wish  that  this  new  imaging  method  produce  an  output  as  in, 
say,  Figure  4.11,  where  the  numerical  value  at  each  point  in  ihe  reconstructed 
field  represents  the  degree  to  which  there  seems  to  be  a  scatterer  at  that  point. 

The  obvious  approach  seems  to  be  to  somewhat  repeat  the  concept  of  the 
pattern-match  method.  That  is,  we  first  assume  that  there  is  a  single  fictitious 
scatterer  (actually,  a  scattering  center)  at  some  specific  location  whose  scattering 
pattern  exhibits  the  frequency  dependence  given  by  Eq.  (5.18).  Second,  we  test 
this  hypothesis  by  using  Eq.  (5. IS)  to  shift  the  information  actually  received, 
from  one  frequency  to  another.  If  the  hypothesis  is  correct,  the  shifted  patterns 
should  line  up  as  seen  in  Figure  5.17;  if  not,  they  should  not  line  up.  And  third, 
the  fictitious  scatterer  is  moved  through  all  locations  in  object  space,  i.e.,  all 
locations  which  could  possibly  contain  an  actual  scatterer.  The  location  at  which 
the  shifted  patterns  became  aligned  the  best  are  then  interpreted  as  the  location 
of  the  true  scatterer. 


159 


The  most  important  issue  which  remains  can  be  expressed  as  a  question: 
what  metric  should  be  used  to  quantify  how  well  the  shifted  patterns  "line  up?” 
Before  this  question  can  be  answered,  we  mus.  first  answer  a  related  one:  do  we 
use  the  complex  received  signal,  or  its  magnitude  (as  was  used  in  Figures  5.15  to 
5.18)? 


There  are  at  least  two  arguments  in  favor  of  avoiding  the  complex  signal. 
First,  for  non-reverberant  scatterers,  the  phase  of  the  received  signals  is  deter¬ 
mined  by  the  wavelength  and  the  path  length.  That  is,  though  the  phase  varies 
with  the  angle  d>R  to  a  receiver,  the  variation  is  due  only  to  the  change  in  the 
path  length  to  the  receiver  —  there  is  no  additional  variation  due  to  some  scat¬ 
tering  phenomenon.  Second,  after  the  received  pattern  at  one  frequency  has  been 
shifted  to  another  frequency,  we  would  need  to  add  some  additional  compensa¬ 
tion  to  the  phase  since  each  point  in  the  shifted  pattern  is  now  at  a  new  location 
in  the  receiver  array  and  therefore  at  a  new  distance  from  the  scatterer.  It  is  not 
clear  at  this  time  if  such  an  additional  compensation  could  be  determined.  For 
these  reasons,  only  the  magnitude  of  the  received  pattern  shall  be  used. 

We  now  return  to  the  original  question  on  determining  how  well  the  shifted 
patterns  line  up.  This  question  is  quite  similar  to  the  one  addressed  by  the 
mismatch  method  developed  in  Section  4.2  if  we  instead  ask  the  question:  how 
well  do  the  shifted  patterns  not  line  up?  To  measure  this,  we  redefine  the  metric 
Zj,  originally  developed  for  comparing  complex  numbers,  to  be 

(7  R. 

i -  kQm))2  (5.i9) 

1  q= 1  i=l  V  7 

where  Nc  is  the  number  of  frequency  pairs  which  are  being  compared.  The 
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two  frequencies  f(q)  and  f(q)  are  any  two  frequencies  at  which  data  has  been 
recorded  (the  appropriate  selection  of  these  frequencies  will  be  investigated  in 
the  next  section).  The  values  /?i  (/(<?))  represent  the  pattern  at  array  locations 
XR  which  have  been  shifted  from  frequency  f(q)  to  frequency  f(q).  The  values 
jRi(7(<7))  also  represent  the  pattern  ai  array  locations  X-1 ;  however,  these  are 
values  which  have  not  been  shifted  in  frequency  but  which  have  instead  remained 
at  their  original  frequency  f(q).  Thus,  if  the  patterns  at  each  pair  of  frequencies 
line  up,  we  have  zj  =  0;  otherwise.  z\  can  grow  to  be  a  very  large  number. 

The  normalization  factor  in  Eq.  (5.19)  is  equal  to  the  number  of  terms 
compared  and  is  ideally  given  by  hrlil  =  Nc  ■  NR  where  NR  is  the  number  of 
receivers.  However,  as  was  seen  in  Figure  5.16,  the  number  of  terms  compared  will 
be  decreased  when  there  are  one  or  more  undetermined  points  after  a  frequency 
shift. 

Each  term  of  Eq.  (5.19)  compares  two  values:  R{  (/(?))  and  R,(^f(q)j.  The 
former  requires  no  calculation  as  it  is  one  of  the  values  directly  measured  at  the  z th 
receiver.  The  latter,  however,  must  be  estimated  as  follows.  First,  we  calculate 
the  shifted  pattern  locations  using  Eq.  (5.18)  with  /  — ►  f(q)  and  f  — ►  f{q).  That 
is,  the  data  values  R{  which  were  originally  recorded  at  receiver  locations  XR  are 
now  considered  to  be  data  values  at  the  shifted  receiver  locations  XR.  Second,  we 
must  use  interpolation  to  evaluated  the  shifted  pattern  at  the  original  receiver 
locations  since  we  only  know  the  values  RiQ{q)^  at  those  locations.  We  can 
write  this  as 

&(/(«))  =  interpolate)  .v/t  n  (A',y  fiq'l'j  for  m=  1,2, ,  Nr 


(5.20) 
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where  the  INTERPOLATE  operation  uses  the  set  of  values  in  its  second  operand  to 
estimate  the  function  at  its  first  operand. 

Once  all  values  of  zj  have  been  i  mputed,  a  final  transformation  is  per¬ 
formed  as  given  in  Eq.  (4.17).  That  is  n  cases  where  the  frequency-shifted 
patterns  line  up  the  best  (i.e.,  there  is  a  r  mum  of  mismatch)  are  mapped  to 
1  while  the  average  case  is  mapped  to  0. 

5.4.3  Choosing  the  Frequencies 

One  of  the  input  parameters  to  the  program  which  generates  the  swept- 
frequency  image  is  a  table  describing  the  pairs  of  frequencies  which  are  to  be 
compared: 

f(q)  and  f(q)  for  q  =  1, 2, . . .  ,  jVC.  (5.21) 

Although  such  a  table  can  accommodate  any  possible  combination  of  frequencies, 
several  pragmatic  considerations  limit  the  combinations.  First,  in  order  to  help 
simplify  the  specification  of  the  table,  we  always  assign  the  higher  frequency  to 
the  second  member  of  the  pair  so  that  f(q)  >  f(q).  Second,  we  always  shift 
the  lower  frequency  to  the  higher  frequency  since  there  will  be  more  sinusoid-like 
cycles  across  the  array  at  the  higher  frequencies  and  thus  a  greater  sensitivity  to 
mismatch. 

To  better  understand  how  the  choice  of  frequencies  affected  the  new  imag¬ 
ing  method,  a  simulation  was  prepared  and  run  for  a  variety  of  frequencies.  All 
of  the  results  reported  in  the  remainder  of  this  section  use  this  same  simulation. 

The  simulation  consists  of  a  source  at  (20",  0)  with  a  pair  of  point  reflectors 
centered  at  (—20”,  100”)  and  separated  by  12"  so  that  the  line  connecting  the 
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reflectors  is  parallel  to  the  X-axis.  An  array  of  64  receivers,  spaced  0.5”  apart, 
is  centered  on  the  X-axis.  Only  a  one-dimensional  reconstruction  along  a  line  at 
Z  —  100"  was  generated  in  these  tests  since  we  were  mainly  interested  in  lateral 
resolution. 

The  first  three  tests,  shown  in  Figures  5.19,  5.20,  and  5.21,  use  a  single  pair 
of  frequencies  arranged  so  that  A /  =  2.5  kHz,  A /  =  5  kHz,  and  A /  =  10  kHz,  in 
that  order.  These  figures  reveal  several  phenomena.  First,  a  rough  measurement 
of  the  distance  between  the  peaks  reveal  that  they  are  separated  by  a  distance 
proportional  to  the  wavelength  of  the  difference  in  frequencies,  c/A/,  rather 
than  the  wavelength  of  either  of  the  two  frequencies  actually  used.  Second,  as 
the  separation  of  the  peaks  becomes  less,  there  are  more  false  peaks.  Third,  the 
peaks  become  more  narrow  as  the  difference  in  frequencies  A /  increases,  even 
thought  the  higher  frequency  (30  kHz)  was  the  same  in  all  cases. 

A  fourth  phenomenon,  which  we  shall  call  a  “missed  comparison,”  arises 
at  either  end  of  the  curve  for  the  largest  A/,  in  Figure  5.21.  This  phenomenon 
was  first  noted  in  Section  5.4.1  and  is  due  to  the  fact  that  the  frequency-shifted 
pattern  has  moved  completely  off  the  array,  making  it  impossible  to  compare 
it  to  the  unshifted  pattern.  So  that  we  can  distinguish  this  phenomenon  from 
that  of  zero  mismatch,  the  imaging  program  assigns  a  value  of  0.01  or  1%  of  the 
maximum  to  the  image  when  this  problem  occurs.  These  values  thus  appear  as 
low-level  plateaus.  These  missed  comparisons  are  more  likely  to  occur  whenever 
(1)  the  value  of  A /  is  large,  and  (2)  the  fictitious  scatter  is  at  a  large  angle 
relative  to  the  Z-axis. 

What  we  would  like  is  to  have  the  high  precision  which  occurs  when  A /  is 
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Figure  5.19.  Swept— frequency  reconstruction  of  2  point  reflectors  at  Z— 100  using  Af— 2.5  kHz  (27.5  and 
30  kHz). 


Swept-frequency  reconstruction  of  2  point  reflectors  at  Z=100"  using  Af=10  kHz  (20  and  30 
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large  with  the  lesser  number  of  false  peaks  which  occurs  when  A /  is  small.  One 
approach  to  this  ideal  would  be  to  use  a  smaller  A/,  but  repeated.  In  Figure 
5.22.  3  pairs  of  frequencies  are  used  where  A /  =  5  kHz  for  each  pair.  The  results 
are  essentially  the  same  as  when  a  single  pair  of  frequencies  was  used  (Figure 
5.20);  thus,  nothing  has  been  gained. 

One  additional  question  which  arose  was  whether  the  frequencies  used 
could  be  lowered  without  sacrificing  image  resolution.  The  case  shown  in  Figure 
5.20  has  been  repeated  for  Figure  5.23,  where  A /  was  kept  at  5  kHz  but  the 
frequencies  have  been  lowered  from  25  kHz  and  30kIIz  to  10  kHz  and  15  kHz. 
While  the  width  of  the  peaks  are  essentially  the  same,  the  number  of  missed 
comparisons  has  grown  considerably  due  to  the  fact  that  A/,  as  a  percentage  of 
the  frequencies,  has  increased. 

The  combinations  of  frequencies  examined  so  far  indicate  that  the  false 
peaks  can  be  reduced  if  a  set  of  values  for  A /  are  chosen  where  each  is  an 
irrational  multiple  of  the  others.  One  convenient  method  for  the  construction  of 
such  a  set  is  to  select  frequencies  with  a  logarithmic  distribution,  where  each  is 
compared  to  the  highest.  Such  frequencies  are  given  by 

fi  =  fmin  +  Umax  ~  fmxn)  *  Iog10  (l  +  (*'  -  for  »'  =  1, 2,  .  .  .  ,  iVF 

(5.22) 

where  the  minimum  frequency  fmin,  maximum  frequency  fma x,  and  the  number 
of  frequencies  NF  have  been  given.  Such  a  set  would  then  be  used  to  make  NF  —  l 
comparisons  since  we  cannot  compare  fmax  to  itself. 


A  set  of  20  frequencies,  from  10  kHz  to  30  kHz,  was  constructed  and  is 
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Figure  5.22.  Swept-frequency  reconstruction  of  2  point  reflectors  at  Z-100"  using  3  pairs  of  frequencies 
having  Af=5  kHz  (15  and  20  kHz,  20  and  25  kHz,  and  25  and  30  kHz). 


scattering  center 


Figure  5.23.  Swept-frequency  reconstruction  of  2  point  reflectors  at  Z=100"  using  Af=5  kHz  at  relatively 
low  frequencies  (10  and  15  kHz). 
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shown  in  Table  5.4.  A  reconstruction  using  this  set,  shown  in  Figure  5.24,  exhibits 
most  of  the  characteristics  we  desire.  Though  the  main  peak  is  no  longer  as  sharp 
as  for  the  A/  =  10  kHz  case,  the  false  peaks  have  been  essentially  eliminated. 
When  the  number  of  frequencies  was  increased  or  decreased,  the  results  (not 
shown  here)  showed  that  the  main  peak  was  unchanged  while  the  small  false  peaks 
decreased  or  increased,  respectively,  by  a  small  amount.  Therefore,  all  of  the 
remaining  swept-frequency  simulations  in  this  chapter  use  the  set  of  frequencies 
shown  in  Table  5.4. 

5.5  Locating  Specular  Scatterers 

The  ability  of  the  swept-frequency  method  to  locate  specular  scatterers 
has  been  explored  through  a  series  of  simulations.  In  every  case,  a  pair  of  point 
reflectors  are  used  which  are  separated  by  12"  along  a  line  parallel  to  the  X-axis. 
Also,  each  simulation  uses  an  array  of  64  receivers  separated  by  0.5"  which  is 
centered  on  the  X-axis,  and  the  20  frequencies  listed  in  Table  5.4. 

The  following  subsections  examine  the  results  when  the  scattering  center 
(the  point  midway  between  the  two  point  reflectors)  is  at  different  ranges;  when 
it  is  at  different  lateral  locations;  when  t  he  scattering  center  and  source  have  no 
offsets  (i.e.,  they  lie  on  the  Z-axis);  and  when  there  are  two  sets  of  scatter  ws. 

All  contour  plots  will  include  lines  at  10f/u,  20%,  . . . ,  90%  of  the  maximum, 
and  will  include  a  sketch  of  the  receiver  array  (a  line  with  a  tic  mark  at  each 
receiver)  and  source  when  possible.  All  hidden-surface  plots  will  have  data  values 
below  0.005  (bad  matches)  erased,  and  values  of  0.01  (missed  comparisons)  shown 
as  plateaus. 
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Table  5.4.  Logarithmic  set  of  frequency  comparisons. 


Low  Frequency  (Hz) 


10.000.00 


13, 36S.09 


15.7Sb.96 


17, 680. OS 


19,232.18 


20,54S.53 


21,691.39 


22.701.21 


23,605.76 


24,424.93 


25,173.46 


25,862.57 


26,501.00 


27,095.71 


27,652.29 


2S, 175.34 


2S,6SS.6S 


29,135.50 


29,578.50 


High  Frequency  (Hz) 


30.000.00 


30,000.00 


30,000.00 


30,000.00 


30,000.00 


30,000.00 


30,000.00 


30,000.00 


30.000.00 


30,000.00 


Difference  (Hz) 


20,000.00 


16,6 31.91 


14,211.04 


12,319.92 


10,767.S2 


9,451.47 


8,308.61 


7,298.79 


6,394.24 


5,575.07 


30,000.00 

4,S26.54 

30,000.00 

4,137.43 

30,000.00 

3,499.00 

30,000.00 

2,904.29 

30,000.00 

2,347.71 

30,000.00 

1,824.66 

30,000.00 

1,311.32 

30,000.00 

S64.50 

30,000.00 

421.50 

Figure  5.24.  Swept-frequency 
distribution  of  20  frequencies 


5.5.1  Specular  Results  at  Different  Ranges 

The  simulations  described  in  this  section  were  selected  to  explore  how  the 
swept -frequency  image  of  a  pair  of  point  reflectors  varies  with  the  range  of  the 
reflectors.  In  all  of  the  simulations  discussed  in  this  section,  the  lateral  location 
of  the  scattering  center  is  X  =  —20". 

The  reconstruction  of  a  pair  of  reflectors  whose  scattering  center  is  at 
(—20". 25")  is  shown  in  Figure  5.25.  An  expanded  view  of  the  results  near  the 
scattering  center,  shown  in  Figure  5.26,  indicates  that  there  is  a  small  error  in 
the  location  of  the  peak.  Such  an  error  might  have  been  expected  since,  at  such 
a  small  range,  the  farfleld  assumptions  built  into  Eq.  (5. IS)  would  no  longer  be 
valid.  Another  view  of  this  reconstruction,  shown  in  Figure  5.27,  exhibits  the 
plateaus  which  were  discussed  earlier  and  which  are  areas  of  missed  comparisons, 
i.e.,  where  the  frequency  shift  of  the  patterns  from  f(q)  to  f(q)  moved  the  patterns 
entirely  off  the  end  of  the  receiver  array. 

Figure  5.27  also  shows  several  sharp  spikes  which  will  be  called  “sparse 
comparisons.”  Whereas  a  missed  comparison  occurs  when  the  number  of  com¬ 
parisons,  A  ,  has  gone  to  zero,  a  sparse  comparison  occurs  when  Nc  has  gone  to 
some  small  but  non-zero  value.  Thus,  even  though  the  patterns  may  not  line  up 
at  all,  there  are  so  few  comparisons  made  that  the  amount  of  mismatch  is  very 
small,  i.e.,  the  amount  of  apparent  match  is  very  large.  This  phenomenon  could 
be  eliminated  by  treating  cases  where  Nc  is  below  some  threshold  as  if  they  were 
cases  of  missed  comparisons;  however,  this  has  not  been  done  in  the  figures  in 
this  chapter  so  that  the  phenomenon  can  be  observed. 


reflector  j  scattering  center 


frequency  reconstruction  of  2  point  reflectors  centered  at  (-20", 25")  with  a  source 
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la  Figure  5. 28,  we  have  doubled  the  range  so  that  the  scattering  center 
is  now  at  (—20", 50").  While  the  error  in  the  location  of  the  peak  has  been 
reduced,  the  spread  of  the  peak  has  increased  —  that  is,  the  resolution  has  been 
decreased.  This  is  even  more  apparent  in  Figure  5.29,  where  we  also  see  that  the 
range  resolution  has  become  much  poorer. 

When  the  range  is  again  doubled  to  Z  —  100",  the  range  resolution  dis¬ 
appears.  Figure  5.30  shows  that  the  lateral  resolution  has  also  become  poorer. 
However,  if  the  lateral  resolution  were  to  remain  the  same  as  the  range  went 
to  infinity  (as  might  be  suggested  by  Figure  5.30),  we  would  have  a  remarkable 
imaging  method.  However,  as  is  shown  in  Figure  5.31  for  a  range  of  Z  =  400", 
this  is  not  the  case  —  the  lateral  resolution  also  decreases  with  range,  as  for  most 
other  imaging  methods.  That  is,  while  the  lateral  resolution  remains  the  same 
as  the  range  of  the  fictitious  source  increases,  it  gets  poorer  as  the  range  of  the 
true  source  increases. 

5.5.2  Specular  Results  at  Different  Bearings 

The  simulations  described  in  this  section  were  selected  to  explore  how  the 
swept-frequency  image  of  a  pair  of  point  reflectors  varies  with  the  lateral  location 
of  the  reflectors.  In  all  of  the  simulations  discussed  in  this  section,  the  range  of 
the  scattering  center  is  Z  =  50". 

We  begin  by  reviewing  Figures  5.2S  and  5.29  from  the  previous  section, 
where  the  lateral  location  of  the  scattering  center  was  X  =  —20"  and  a  certain 
amount  of  range  resolution  was  exhibited.  Next,  we  move  the  scattering  center 
to  X  =  0.  The  results,  shown  in  Figure  5.32,  reveal  that  all  range  resolution  has 
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Figure  5.29.  Swept-frequency  reconstruction  of  2  point  reflectors  centered  at  (-20", 50")  with  a  source 
at  (20", 0),  showing  sparse  comparisons  (spikes)  and  missed  comparisons  (plateaus). 
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Swept-frequency  reconstruction  of  2  point  reflectors  centered  at  (-20’\100")  with  a  source 


Figure  5.31.  Swe 
with  a  source  at 
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1  een  lost,  and  that  the  region  where  a  good  match  was  found  makes  a  strange 
curve  where  Z  <  3'J".  Figure  5.33  shows  that  the  missed  and  sparse  comparisons 
are  also  still  present,  and  that  the  match  has  a  very  minor  maximum  at  the  true 
scattering  center. 

When  the  scattering  center  is  again  moved  laterally,  to  A’  =  'JO".  Figures 
5.34  and  5.35  sh-wn  that  the  results  are  quite  similar  to  the  complementary  case 
where  tiie  scattering  center  was  at  A  =  —20"  (see  Figures  5. 28  and  5.20). 

5.5.3  Specular  Results  with  No  Offsets 

A  single  simulation  is  used  in  this  section  to  examine  how  the  the  swept- 
frequency  method  reacts  to  the  case  where  there  are  no  offsets  in  the  location 
of  the  source  and  scattering  center  —  that  is,  they  all  lie  on  the  normal  to  the 
receiving  array.  The  results  of  this  simulation,  shown  in  Figures  5.36  and  5.37, 
indicate  that  the  swept-frequencv  method  has  great  difficulty  in  determining  the 
range  of  the  scattering  center  in  this  situation. 

The  result  shown  in  Figures  5.36  and  5.37  might  have  been  expected  since 
the  frequency-dependent  shift  of  the  patterns  can  be  shown  to  be  about  the 
same,  as  the  range  of  the  fictitious  source  changes.  What  was  not  expected  was 
the  equally  poor  result,  seen  in  the  previous  section  in  Figures  5.32  and  5.33, 
when  the  source  was  to  the  side,  at  Z  =  20". 

5.5.4  Results  Using  Two  Sets  of  Specular  Scatterers 

A  single  simulation  is  used  in  this  section  to  examine  how  the  swept- 
frequency  method  reacts  to  the  case  where  there  are  two  pairs  of  reflectors.  The 


Figure  5.33.  Swept-frequency  reconstruction  of  2  point  r< 
(20'\0),  showing  sparse  comparisons  (spikes)  and  missed  c 
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Figure  5.35.  Swept-frequency  reconstruction  of  2  point  reflectors  centered  at  (20",50  )  with  a  source  at 
(20", 0),  showing  sparse  comparisons  (spikes)  and  missed  comparisons  (plateaus). 


Figure  5.36.  Swept-f  requency  reconstruction  of  2  point  reflectors  centered  at  (0,50'')  with  a  source  a 
(0,0). 


results  shown  in  Figure  5.3S  indicate  where  the  reflectors  have  been  located.  In 

this  figure,  we  have  shown  only  that  portion  of  the  reconstructed  field  near  the 

/ 

reflectors  and  have,  in  addition,  renormalized  the  data  to  the  maximum  within 
this  portion. 

The  entire  reconstructed  field  is  shown  in  Figure  5.39.  Here,  the  size  of  the 
sparse-comparison  spikes  is  so  much  greater  than  that  of  the  peaks  around  the 
scattering  centers  that  the  latter  have  been  made  quite  small  by  normalization. 

A  hidden-surface  plot  of  that  portion  of  the  reconstructed  field  near  the 
reflectors  is  shown  in  Figure  5.40,  with  the  data  again  renormalized  to  the  max¬ 
imum  within  the  portion  shown.  This  figure  indicates  that  the  4  reflectors  have 
formed  multiple  scattering  centers,  among  the  various  pairs  of  reflectors,  with  the 
peak  near  X  =  —10"  perhaps  representing  the  sum  of  2  scattering  centers  where 
one  is  due  to  the  inner  pair  of  reflectors  and  one  is  due  to  the  outer  pair. 

5.6  Locating  Extended  Scatterers 

The  simulations  performed  in  Section  5.5  are  now  exactly  repeated  except 
that  the  two  point  reflectors  are  replaced  with  a  plate  reflector  which  is  12" 
wide  and  which  is  simulated  using  a  summation  of  200  terms.  The  logarithmic 
distribution  of  frequencies  given  in  Table  5.4  is  used  in  each  case. 

5.6.1  Plate  Results  at  Different  Ranges 

The  simulations  described  in  this  section  were  selected  to  explore  how  the 
swept- frequency  image  of  a  plate  reflector  varies  with  the  range  of  the  reflector. 
In  all  of  the  simulations  discussed  in  this  section,  the  lateral  location  of  the  center 
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of  the  plate  is  A"  =  —20". 

The  reconstruction  of  a  plate  whose  center  is  at  (—20", 25")  is  shown  in 
Figure  5.41  and  exhibits  the  same  small  error  in  the  location  of  the  peak  as  was 
seen  in  the  case  of  the  point  reflectors.  It  is  interesting  to  also  note  that  the  range 
ambiguity  extends  along  a  line  which  is  parallel  to  the  Z-axis,  rather  than  along  a 
line  running  radially  from  the  center  of  the  receiver  array.  A  hidden-surface  plot 
of  this  same  information,  Figure  5.42,  has  a  region  showing  a  good  match  along 
the  line  at  Z  —  1",  between  X  =  0  and  X  =  17".  It  is  not  clear  whether  these  are 
merely  a  large  number  of  closely  spaced  sparse-comparison  spikes  or  some  other 
phenomenon.  For  example,  it  is  shown  in  the  next  section  that  such  regions  can 
occur  when  the  array  is  within  the  actual  source’s  central  lobe  —  an  equivalent 
phenomenon  may  be  occurring  with  the  fictitious  source  in  Figure  5.42  . 

In  Figures  5.43  and  5.44,  the  range  has  been  doubled  to  Z  =  50";  as 
expected,  the  range  resolution  has  become  poorer  while  the  small  error  in  the 
main  peak  has  been  eliminated.  When  the  range  is  again  doubled,  to  Z  —  100", 
all  range  resolution  is  lost  as  is  seen  in  Figure  5.45. 

5.6.2  Plate  Results  at  Different  Bearings 

The  simulations  described  in  this  section  were  selected  to  explore  how  the 
swept-frequency  image  of  a  plate  varies  with  the  lateral  location  of  the  plate.  In 
all  of  the  simulations  discussed  in  this  section,  the  range  of  the  plate  is  Z  —  50". 

We  can  consider  the  series  discussed  in  this  section  as  beginning  with 
Figures  5.43  and  5.44,  of  the  previous  section,  as  they  are  for  a  plate  centered  at 
(—20", 50").  When  the  plate  is  then  moved  laterally,  to  A'  =  0,  Figures  5.46  and 


Figure  5.41.  Swept-frequency  reconstruction  of  a  12"  plate  centered  at  (-20", 25")  with  a  source  at 
(20", 0). 
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Figure  5.42.  Swept- frequency  reconstruction  of  a  12  plate  centered  at  (  20  ,25  )  with  a  source  at 
(20”, 0),  showing  sparse  comparisons  (spikes)  and  missed  comparisons  (plateaus). 
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Figure  5.46.  Swept-frequency  reconstruction  of  a  12"  plate  c 
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•5.47  show  that  in  addition  to  a  loss  of  all  range  resolution  (as  was  experienced  for 
the  point  reflectors),  a  major  loss  of  lateral  resolution  is  also  experienced.  When 
the  center  of  the  plate  is  moved  still  further  to  the  right,  to  A'  =  20",  Figures 
5.4S  and  5.49  show  that  we  do  not  see  a  recovery  of  range  or  lateral  resolution  as 
was  seen  tor  the  point  reflectors  —  the  results  continue  to  be  extremely  poor. 

We  can  gain  some  insight  as  to  the  cause  of  these  poor  results  for  this  last 
case,  where  the  center  of  the  plate  was  at  (20",  50"),  by  examining  the  pattern 
recorded  at  the  receiving  array  for  two  frequencies.  These  patterns,  shown  in 
Figure  5.50,  reveal  that  we  are  receiving  that  portion  of  the  plate's  scattering 
pattern  near  the  central  lobe,  where  the  amplitude  varies  rapidly.  Therefore, 
when  Eq.  (5.18)  is  applied  to  shift  the  lower  frequency  to  the  higher  frequency, 
a  poor  match  results  even  when  the  fictitious  scatterer  is  at  the  location  of  the 
true  scatterer.  The  results  when  the  center  of  the  plate  was  further  away  from 
the  source,  as  shown  in  Figures  5.43  and  5.44,  were  much  better  since  we  were 
then  receiving  that  portion  of  the  the  plate’s  scattering  pattern  which  was  far 
away  from  the  central  lobe,  where  the  amplitude  was  not  changing  as  drastically. 

5.6.3  Plate  Results  with  No  Offsets 

As  might  have  been  expected  after  the  previous  section,  when  the  plate 
and  source  are  both  on  the  Z-axis,  as  is  shown  in  Figures  5.51  and  5.52,  the  image 
exhibits  very  poor  resolution.  In  fact,  Figure  5.52  indicates  that  the  match  for 
a  fictitious  scatterer  at  the  location  of  the  true  scatterer  is  actually  a  bit  worse 
(the  trough)  than  for  locations  on  either  side. 


center 


Figure  5.47.  Swept-frequency  reconstruction  of  a  12"  plate  centered  at  (0,50")  with  a  source  at  (20", 0), 
showing  sparse  comparisons  (spikes)  and  missed  comparisons  (plateaus). 
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Figure  5.49.  Swept— frequency  reconstruction  of  a  12"  plate  centered  at  (20  ,50  )  with  a  source  at 
(20", 0),  showing  missed  comparisons  (plateaus). 
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Figure  5.50.  Patterns  at  an  array  on  the  X-axis  (64  receivers)  for  30  kHz  (solid  line)  and  20  kHz 
(dashed  line),  due  to  a  12”  plate  centered  at  (20”, 50”)  and  a  source  at  (20”, 0). 


Figure  5.51.  Swept— frequency  reconstruction  of  a  12"  plate  centered  at  (0,50")  with  a  source  at  (0,0). 
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5.6.4  Results  Using  Two  Plate  Scatterers 

The  results  of  a  simulation  with  two  12"  plates,  centered  at  (  —  20", 50") 
and  (0,50"),  are  shown  in  Figures  5.53  and  5.54.  While  earlier  simulations  might 
have  suggested  that  we  would  not  be  able  to  resolve  the  plate  on  the  Z-axis,  it  is 
not  clear  why  the  offset  plate  cannot  be  located.  It  is  possible  that  the  off-axis 
contributions  of  the  offset  plate  were  essentially  lost  since  they  were  so  much 
lower  than  the  central  lobe  of  the  on-axis  plate. 

5.7  Comparison  to  Beam  Forming 

In  addition  to  the  simulations  of  the  previous  sections,  the  swept-frequency 
method  can  be  evaluated  by  comparing  it  to  the  results  obtained  using  a  well- 
known  method  such  as  beamforming.  The  swept-frequency  result  for  two  point 
reflectors,  12"  apart  and  centered  at  (—20",  1000"),  is  shown  as  the  solid  curve 
in  Figure  5.55.  The  source  is  at  (20", 0).  When  the  signals  from  the  same  two 
reflectors  were  analyzed  by  a  beamformer,  the  curve  in  Figure  5.55  with  the 
short  dashes  resulted.  While  this  curve  has  a  slightly  sharper  peak  than  the 
swept-frequency  method,  it  is  also  in  error  by  several  inches. 

Part  of  the  reason  for  the  error  in  the  beamformer  result  is  probably  due 
to  the  fact  that  it  assumes  that  it  is  receiving  plane  waves.  The  assumption 
would  be  slightly  more  valid  if  one,  instead  of  two,  reflectors  were  used.  When 
this  case  was  simulated,  the  curve  in  Figure  5.55  with  the  long  dashes  resulted, 
and  exhibited  a  peak  which  was  at  the  correct  location  and  slightly  sharper  than 
the  curve  from  the  swept-frequency  method.  A  swept-frequency  curve  using  only 
one  reflector,  for  comparison,  would  be  useless  because  a  single  reflector  exhibits 


Figure  5.53.  Swept-frequency  reconstruction  of  two  12"  plate 
source  at  (0,0). 


reflector 


no  frequency  dependence  whatsoever. 


The  range  of  the  reflectors  used  for  the  previous  figure  was  made  relatively 
large  so  as  to  approach  the  plane-wave  requirement  of  the  beamformer.  When 
the  range  is  reduced  to  100",  the  results,  shown  in  Figure  5.56,  show  that  while 
the  peak  from  the  swept-frequency  method  has  become  much  more  precise,  the 
beamformer  results  have  become  useless. 


5.8  Results  With  Range  Resolution 


In  most  imaging  systems  which  use  multiple  frequencies,  the  received  sig¬ 
nal  is  correlated  with  the  transmitted  one  to  determine  the  time-of-flight  and 
therefore  the  range  of  the  scatterer.  Ermert  and  Karg28  showed  that  the  ambi¬ 
guity  function  for  differences  in  range  A*  near  a  scatterer  at  range  2  was  given 
bv 


A'(A2) 


sin  (A k  •  Az  ■  -fi) 
Ale  ■  Az  ■ 


(5.23) 


where  R  is  the  distance  between  the  center  of  the  receiving  array  and  the  scatterer. 


It  would  be  interesting  to  roughly  simulate  how  such  a  correlation  could 
improve  a  swept-frequency  image,  using  the  same  frequencies  which  are  already 
being  used  in  the  frequency  sweep.  This  can  be  done  rather  trivially  by  simply 
multiplying  any  of  the  normalized  reconstructions  shown  previously  in  this  chap¬ 
ter  by  the  values  provided  by  Eq.  (5.23).  In  particular,  the  reconstruction  shown 
in  Figure  5.30  for  a  pair  of  reflectors  centered  at  (  —  20",  100”)  would  be  a  good 
candidate  as  it  exhibits  good  lateral  resolution  but  no  range  resolution. 

Before  proceeding,  the  definition  of  “bandwidth”  needs  to  be  addressed. 


In  Ermert  and  Karg,  a  continuous  sweep  of  frequencies  were  used  where  the 
bandwidth  A /  was  defined  as  the  difference  between  the  maximum  and  minimum 
frequencies,  with  A k  =  2tt  A  f/c.  However,  in  Figure  5.30.  a  continuous  sweep 
is  not  used,  and  the  difference  between  the  frequencies  being  compared  can  vary 
from  a  low  of  421.5  IIz  to  a  high  of  20,000  Hz.  Therefore,  each  of  these  frequency 
differences  were  used  for  A /  and  the  results  for  the  low  and  high  differences 
shown  in  Figures  5.57  and  5.5S,  respectively. 

5.9  Swept-Frequency  Execution  Times 

Most  of  the  swept-frequency  processing  in  this  chapter  evaluated  the  re¬ 
construction  field  at  101  points  along  the  X-axis  and  at  100  points  along  the 
Z-axis,  for  a  total  of  10,100  points.  In  all  cases,  20  frequencies  and  19  frequency 
comparisons  were  used. 

The  MicroVAX  II  processor  required  about  11,200  seconds  to  evaluate  these 
images  when  two  point  reflectors  were  used;  and  about  11,400  seconds  when  the 
plate  reflector  was  used.  It  is  not  clear  why  the  time  to  process  the  signals  from 
the  plate  required  slightly  more  time  since  the  time  to  simulate  the  scattering  is 
not  included  in  the  elapsed  time. 

These  times  indicate  that  it  took  approximately  1.109  or  1.129  seconds  to 
reconstruct  each  field  point  when  the  point  reflectors  or  plate  reflector  was  used, 
respectively.  These  times  are  on  the  order  of  3  times  longer  that  that  needed  by 
the  pattern-match  or  mismatch  method  to  form  an  image  point  with  a  barrier 
present. 


Figure  5.57.  Swept-frequency  reconstruction  of  2  point  reflectors  centered  at  (-20",  100"),  for  a  source 
at  (20”,0),  with  simulated  range  resolution  for  a  bandwidth  of  421.5  Hz. 
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Figure  5.58.  Swept-frequency  reconstruction  of  2  point  reflectors  centered  at  (-20”,  100"),  for  a  source 
at  (20”, 0),  with  simulated  range  resolution  for  a  bandwidth  of  20,000  Hz. 
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Chapter  6 

Summary  and  Conclusions 

A  study  of  the  major  approaches  to  the  generation  of  acoustical  images 
has  led  to  the  development  of  several  new  methods,  and  some  new  observations 
concerning  the  localization  of  acoustic  scatterers.  The  following  sections  will 
summarize  these  results  and  suggest  further  lines  of  investigation. 

6.1  Basic  Issues 

Three  basic  imaging  methods  have  been  studied:  holography,  beamform¬ 
ing,  and  receiver  correlations.  Of  these,  the  holography  method  seems  to  have 
the  greatest  general  usefulness  since  a  beamformer  will  fail  when  the  wavefronts 
are  not  planar  and  the  correlation  method  will  fail  when  there  is  more  than  one 
source  or  scatterer.  Another,  lesser  known  imaging  method,  built  upon  the  con¬ 
cept  of  a  fictitious  scatterer,  was  investigated  and  shown  to  have  several  powerful 
capabilities  and  the  ability  to  incorporate,  as  special  cases,  the  three  previously 
studied  methods. 

It  was  found  that  the  proper  distribution  of  receiver  elements  in  sparse 
arrays  can  definitely  improve  the  array’s  performance.  However,  there  does  not 
seem  to  be  any  further  improvements  to  be  had  when  the  array  is  already  fully 
populated,  with  receivers  every  half-wavelength. 

A  study  of  the  use  of  complex  signals  to  represent  the  information  at  the 
receivers  revealed  that  significant  range  information  is  lost  with  such  a  represen¬ 
tation,  but  could  be  recovered  by  recasting  an  imaging  method  to  instead  use  the 
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time-of-  flight. 

6.2  High  Resolution  Methods 

Two  categories  of  high  resolution  array  processing  methods  were  exam¬ 
ined:  the  autoregressive-moving  average  (ARMA)  methods;  and  the  minimum 
energy  (ME)  method  and  its  variations  including  linear  predictors  (LP)  and  eigen¬ 
vector  analysis.  While  each  holds  some  promise  of  enhancing  the  resolution,  each 
also  has  some  major  obstacles  which  need  to  be  investigated  further  before  the 
method  can  be  applied  to  the  location  of  scatterers. 

It  does  not  seem  possible  to  enhance  either  category  of  high  resolution 
imaging  to  incorporate  a  priori  information  about  the  scatterers  or  their  relative 
location.  Nor  does  it  seem  possible  to  use  the  full  time-history  of  the  received 
signals  so  as  to  avoid  the  phase  ambiguity  described  in  Section  2.6.  On  the 
other  hand,  the  two  imaging  methods  investigated  with  simulations  in  this  study 
(pattern-match  and  swept-frequency)  could  be  modified  to  use  such  information. 
It  would  be  interesting  to  investigate  whether  such  a  modification  would  yield  a 
higher  resolution  than  ARMA  or  ME  methods. 

6.2.1  ARMA  Methods 

The  major  question  in  applying  this  method  is  whether  any  extrapolation 
can  predict  a  scattering  function  at  angles  beyond  those  subtended  by  the  re¬ 
ceiving  array.  While  this  is  probably  possible  for  simple  scatterers,  further  work 
needs  to  be  done  to  determine  if  robust  extrapolation  is  possible  for  the  general 
case.  In  the  absence  of  prior  knowledge  about  the  scatterer,  it  seems  that  these 
methods  should  be  used  with  great  caution.  A  review  of  the  literature  suggests 
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that  this  sort  of  problem  has  not  been  considered  since  most  workers  have  been 
seeking  only  the  location  of  simple  plane-wave  sources. 

In  those  cases  where  an  ARM  A  spectral  estimator  can  properly  extrapo 
late  data  beyond  the  receiving  array,  the  estimator  would  provide  an  enhanced 
resolution  even  in  the  absence  of  noise. 

6.2.2  ME  Methods 

The  major  obstacle  in  applying  minimum  energy  methods  is  their  inability 
to  work  with  any  wavefront  shape  other  than  ones  that  are  planar.  Thus,  while 
intuition  would  predict  that  an  imaging  method  should  yield  better  results  as 
the  object  being  imaged  moved  closer,  these  methods  instead  progress  towards 
uselessness  as  the  object  moves  closer.  While  both  methods  can  be  modified 
to  accommodate  any  specific  wavefront,  no  adaptive  method  is  known  for  the 
determination  of  what  this  wavefront  should  be  —  it  must  be  chosen  by  some 
other  means.  While  this  problem  does  not,  technically  speaking,  rule  out  these 
methods,  it  means  that  they  can  be  used  only  to  locate  scatterers  which  are  in 
the  farfield.  No  scatterer,  no  matter  how  simple  its  scattering  function,  can  be 
handled  if  it  is  too  close. 

A  great  deal  of  progress  has  been  made  in  accommodating  coherent  sources 
during  the  last  few  years  and  so  this  may  no  longer  be  an  obstacle  in  applying  the 
ME  method  to  scatterers.  The  observation  by  Duckworth41  that  small,  naturally 
occurring  motions  can  provide  the  needed  decorrelation  is,  on  the  one  hand,  a 
promising  solution  although  it  has  yet  to  be  established  whether,  on  the  other 
hand,  such  motion  will  smear  the  resulting  image  to  the  extent  that  the  high 
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resolution  method’s  enhancements  are  lost. 

6.3  Pattern-Match  and  Mismatch  Imaging- 

Two  new  imaging  methods,  called  the  pattern-match  method  and  the  mis¬ 
match  method,  were  developed  from  the  concept  of  a  fictitious  scatterer.  The 
mismatch  method  is  similar  to  the  pattern-match  method  but  uses  a  new  ap¬ 
proach  to  the  comparison  of  two  complex  signals.  The  following  sections  review 
the  work  which  has  been  done  in  this  study,  and  then  suggests  several  questions 
and  variations  which  may  make  these  new  methods  yet  more  useful. 

6.3.1  Review  of  the  Simulations 

It  has  been  shown  that,  to  the  extent  that  the  simulation  of  the  diffrac¬ 
tion  phenomena  is  valid,  the  new  methods  are  powerful  tools  for  the  analysis  of 
pressure  waves  which  have  encountered  a  known  disturbance.  These  methods 
can  provide  a  superior  reconstruction  to  the  holography  method  in  many  cases 
when  free-space  conditions  do  not  exist.  The  application  to  experimental  data  is 
limited  only  by  the  accuracy  of  the  simulation,  of  the  experimental  diffraction  or 
other  propagation  phenomena,  which  is  used  in  the  estimate  of  the  signal  from 
the  fictitious  source. 

To  test  the  new  methods,  a  computer  simulation  of  diffraction  by  a  barrier 
was  implemented.  This  simulation  automatically  analyzed  the  relative  locations 
of  the  source,  barrier,  and  receiver,  and  correctly  determined  in  which  of  the 
3  possible  regions  the  receiver  was  located.  Though  developed  for  sources  and 
receivers  far  from  the  barrier,  the  simulation  was  sufficiently  accurate  for  the 
holography  method  to  form  the  image  of  a  reflection  when  the  source  was  in 
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front  of  the  barrier. 

The  simulations  revealed  that,  when  the  acoustic  source  is  partially  hidden 
from  the  receiving  array  by  a  barrier,  the  pattern-match  method  reports  that  the 
source  is  at  some  location  behind  the  barrier,  without  reporting  an  exact  location; 
on  the  other  hand,  the  holography  method  falsely  reports  that  the  source  is  at 
the  edge  of  the  barrier. 

The  mismatch  method  excludes  more  locations,  from  being  possible  source 
positions,  than  the  pattern-match  or  holography  methods.  However,  it  is  not 
clear  that  the  resulting  spike-like  response  would  always  be  more  desirable  — 
without  the  pattern-match  image  as  a  guide,  these  spikes  could  be  interpreted  as 
additional  sources. 

It  was  shown  that  the  new  methods  are  not  significantly  affected  by  small 
errors  in  their  knowledge  of  the  environment,  such  as  an  error  in  the  range  of  a 
barrier. 


Although  it  is  not  clear  whether  the  new  methods  can  be  considered  to 
be  linear  processors,  simulations  using  two  sources  indicate  that  they  at  least 
approximate  the  linear  response  of  the  holography  method. 

The  new  methods  are  generally  slower  than  the  holography  method;  how¬ 
ever,  it  seems  that  this  should  not  be  considered  to  be  an  obstacle  at  this  stage 
as  ever-faster  computer  system  may  make  this  issue  irrelevant. 


6.3.2  Future  Work:  Composite  Sources 


In  this  study,  only  the  pattern  from  a  single  fictitious  source  has  been 
used  in  the  search  for  the  true  source,  even  in  cases  when  there  was  in  fact  two 
sources.  However,  when  a  priori  information  can  establish  the  number  of  sources 
and  their  relative  orientation,  we  could  treat  this  set  of  sources  as  one  composite 
fictitious  source,  and  estimate  the  pattern  which  would  have  been  received  at  the 
array  due  to  each  possible  location  of  the  composite  fictitious  source.  This  might 
be  extremely  effective  when  the  object  being  sought  can  be  characterized  by  a 
set  of  specular  reflections  whose  relative  separation  and  strength  is  known. 

This  same  approach  could  be  taken  further  to  include  the  case  of  an  ex¬ 
tended  source  or  scatterer.  We  would  continue  to  need  an  algorithm  which  can 
estimate  the  pattern  at  the  receiving  array  for  each  possible  location  of  the  ex¬ 
tended  fictitious  object.  This  algorithm  could  either  be  purely  theoretical,  or 
could  perform  the  appropriate  interpolation  upon  a  sufficiently  extensive  set  of 
empirical  measurements. 

6.3.3  Future  Work:  Noise  and  Interference 

Although  this  study  has  laid  the  groundwork  for  the  use  of  the  pattern- 
match  and  mismatch  methods,  actual  applications,  such  as  the  guidance  of  a 
mobile  robot,  will  doubtlessly  encounter  a  variety  of  sources  of  noise  and  inter¬ 
ference. 

An  investigation  of  these  problems  should  begin  by  introducing  isotropic 
noise  into  the  simulations;  if  this  reveals  a  problem,  correlations  between  elements 
might  be  used  in  place  of  the  simple  phase  calculations.  In  addition,  directional 
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noise  (representing,  say,  a  fan  in  the  robot’s  vicinity)  should  also  be  investigated. 

Other  environmental  problems  which  might  be  considered  could  include 
additional  barriers,  errors  in  the  knowledge  of  the  receiver  array  elements,  and 
the  occurrence  of  motion  during  signal  acquisition. 

6.3.4  Future  Work:  Mismatch  Variations 

There  would  seem  to  be  an  almost  unlimited  number  of  variations  of  the 
mismatch  method. 

In  this  study,  the  mismatch  was  equated  to  the  “distance”  between  two 
complex  numbers;  this  mismatch  metric  thus,  more  or  less,  gives  the  same  weight 
to  differences  in  phase  as  it  does  to  differences  in  amplitude.  These  weights 
could,  for  example,  be  varied  by  calculating  the  phase  and  amplitude  differences 
separately  and  then  weighing  them  in  some  other  manner.  In  some  cases,  it  may 
be  useful  to  allow  the  sign  of  the  differences  to  have  an  effect. 

In  any  of  these  variations,  the  resulting  match  could  also  be  raised  to 
some  higher  power,  or  inverted,  so  as  to  give  a  greater  emphasis  to  the  large 
mismatches. 

In  addition,  any  of  these  variations  could  be  normalized  to  the  amplitude 
of  one  of  the  complex  numbers  so  that,  for  example,  mismatches  between  small 
numbers  are  weighed  more  heavily  than  equal-sized  mismatches  between  large 
numbers. 

In  would  quite  useful  to  investigate  how  the  mismatch  method  achieves 
its  improvements,  in  those  regions  where  it  offers  improvements  over  the  pattern- 
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match  method  or  holography,  with  the  goal  of  extending  its  improvements  to  all 
cases.  In  particular,  it  may  be  that  a  diversity  of  frequencies  could  eliminate 
the  false  peaks,  such  as  does  a  diversity  of  values  for  A/  in  the  swept-frequency 
method. 

6.3.5  Future  Work:  Time-of-Flight 

In  this  study,  both  match  methods  represent  the  time  delay  seen  at  each 
receiver  as  a  complex  number.  As  discussed  in  Section  2.6,  when  the  transmitted 
signal  is  broadband,  considerably  more  information  is  present  than  can  be  rep¬ 
resented  by  a  single  complex  number.  In  these  cases,  it  is  possible  to  adapt  the 
pattern-match  methods  to  retain  and  use  the  full  time-of- flight  (TOF)  informa¬ 
tion  content  of  the  received  signal. 

To  review,  in  the  method  used  in  Chapter  4,  the  complex  signal  from 
the  single-frequency  fictitious  source  was  estimated  at  each  receiver  location. 
However,  it  would  be  just  as  easy  to  express  this  as  a  complex  amplitude  plus  a 
TOF.  The  path  length  and  hence  the  TOF  is  already  present  in  the  terms  like 
exp(y  k  . . .)  in  the  expressions  for  pj,  pi,  and  pr .  This  allows  us  to  calculate  the 
time-history  to  be  expected  at  each  receiver  location  from  the  fictitious  source.  A 
comparison  of  the  actual  signal’s  time-history  with  that  of  the  fictitious  source, 
summed  across  the  receiver  array,  would  then  yield  the  match  for  that  location 
of  the  fictitious  source.  An  investigation  of  such  an  imaging  method  would  also 
need  to  consider  techniques  for  estimating  the  transmitter’s  waveform  and  time 
of  transmission  when  those  quantities  are  not  known. 
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6.4  Swept-Frequency  Imaging 

A  new  imaging  method  which  uses  multiple  frequencies  was  developed  and 
tested  using  simulated  reflections  from  point  reflectors  and  a  plate.  The  following 
sections  review  the  new  method,  the  results  of  the  simulations,  and  a  variety  of 
developments  which  may  make  the  new  method  yet  more  useful. 

6.4.1  Review  of  Swept-Frequency  Imaging 

A  new  imaging  technique  has  been  developed  which  for  the  first  time  ex¬ 
ploits  the  fact  that  an  object's  scattering  pattern  has  a  deterministic  dependence 
upon  the  frequency  of  the  acoustic  wave  which  it  is  scattering.  This  development 
also  borrowed  the  central  concept  of  the  mismatch  method  to  form  a  metric  which 
indicates  whether  a  given  scattering  pattern  exhibits  the  expected  frequency  de¬ 
pendence. 

To  test  the  new  method,  computer  simulations  of  two  categories  of  scat¬ 
tered  were  implemented  and  tested:  specular  scattered,  represented  by  point 
reflectors:  and  extended  scattered,  represented  by  a  reflective  plate.  The  latter, 
though  implemented  as  a  rather  simple  summation,  proved  to  be  highly  accurate, 
yielding  an  unexpected  but  entirely  appropriate  representation  of  the  scattering 
pattern  at  high  frequencies  and  large  plate  dimensions. 

The  simulations  reveal  that  the  swept-frequency  method  can  provide  good 
range  and  lateral  resolution  in  many  circumstances,  but  can  loose  resolution  when 
the  scatterer  is  moved  to  the  normal  to  the  array  of  receivers,  or  to  larger  ranges. 
The  reconstruction  of  a  plate  can  also  deteriorate  when  it  is  moved  so  that  the 
array  of  receivers  is  within  that  portion  of  the  plate’s  scattering  pattern  where 


the  amplitude  varies  rapidly. 


When  multiple  scatterers  were  present,  the  simulations  showed  that  the 
swept-frequencv  method  can  yield  poor  results.  Multiple  reflectors  caused  the 
image  to  have  a  peak  at  the  scattering  center  formed  by  every  possible  pair  of 
reflectors.  Multiple  plates  yielded  a  rather  useless  image,  though  the  results  may 
have  been  partially  due  to  the  rapidly  varying  amplitude  noted  previously. 

It  was  shown  that  a  simple  range  calculation,  as  is  often  done  with  broad¬ 
band  pulses,  could  be  combined  with  the  swept-frequency  method  to  provide  an 
enhanced  range  resolution. 

A  comparison  with  a  beamformer  was  performed  which  indicates  that, 
while  the  beamformer  may  provide  a  higher  resolution  at  very  large  ranges  (so 
that  it  sees  only  plane  waves),  the  swept-frequency  method  provided  comparable 
or  superior  resolution  at  all  ranges. 

6.4.2  Future  Work:  Amplitude  Dependence 

The  poor  results  experienced  when  attempting  to  image  tiie  plate  at  cer¬ 
tain  locations  was  shown  to  be  due  to  the  plate’s  scattering  pattern  changing 
amplitude,  as  well  as  angle,  as  the  frequency  was  varied.  A  similar  problem  could 
occur  with  other  extended  scatterers. 

It  may  be  possible  to  avoid  this  problem  if  we  replaced  the  pattern  across 
the  array  at  each  frequency  with  an  artificial  pattern  designed  so  that  its  peaks 
and  nulls  are  at  the  same  location  as  the  original.  In  other  words,  all  patterns 
would  be  modified  to  be  similar  to  those  obtained  from  point  reflectors.  It  may 
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even  be  possible  to  skip  the  actual  generation  of  the  artificial  pattern  and  use 
only  a  table  of  peak  and  null  locations  since  the  interpolation,  which  must  be 
done  anyway,  will  be  implicitly  fitting  such  a  pattern  to  its  input  table  of  values. 

One  of  the  challenges  of  such  a  development  would  be  the  creation  of  an 
algorithm  which  could  find  the  desired  locations  in  the  scattering  pattern  of  a 
plate  near  its  central  lobe  —  that  is,  the  phenomena  we  are  looking  for  are  merely 
local  peaks  and  nulls. 

6.4.3  Future  Work:  Phase  Dependence 

It  is  sometimes  possible  to  treat  the  location  of  the  scatterer  as  the  source 
of  a  new  spherical  wave.  Then,  since  we  know  the  wavelength  at  each  frequency, 
the  location  of  the  receivers,  the  location  of  the  transmitter,  and  the  location  of 
the  fictitious  scatterer,  it  would  be  possible  to  determine  the  relative  phase  which 
we  would  expect  to  see  at  each  receiver.  A  more  sensitive  measure  of  mismatch 
could  thus  be  constructed  by  also  comparing  the  expected  and  actual  phases  at 
each  receiver  and  frequency.  Such  an  algorithm  would,  in  effect,  be  looking  for 
spherical  waves  in  the  received  information  which  correspond  to  the  location  of 
the  fictitious  scatterer. 

An  investigation  into  this  enhancement  would  need,  however,  to  first  con¬ 
sider  whether  the  phase  of  the  scattered  wavefront  acts  like  that  of  a  spherical 
wave  centered  at  the  scatterer  —  that  is,  does  the  scattering  process  result  in  a 
scattering  phase  pattern  as  well  as  a  scattering  amplitude  pattern? 

A  yet  more  ambitious  variation  of  this  approach  could  also  be  investigated 
where  the  estimate  of  the  phase  was  replaced  with  an  estimate  of  tire  time-of- 


(light.  This  would  seem  to  also  require  that  the  simple  comparison  of  phases  be 
replaced  with  some  manner  of  correlation. 

6.4.4  Future  Work:  Noise  and  Interference 

In  Section  6.3.3,  it  was  noted  how  various  environmental  events  could 
cause  a  corruption  of  the  pattern-match  and  mismatch  methods.  These  same 
considerations  apply  to  the  swept-frequencv  method  in  some  applications,  such 
as  for  the  guidance  of  a  mobile  robot. 

6.4.5  Future  Work:  Use  of  Frequencies 

There  are  several  questions  concerning  the  use  of  frequencies  which  might 
be  investigated.  It  would  be  useful  to  determine  an  analytical  or  other  approach 
to  the  selection  of  frequencies  which  yielded  the  maximum  resolution  while  also 
eliminating  false  peaks.  These  goals  may  also  be  better  met  by  weighing  the 
results  of  the  frequency  comparisons  unequally,  instead  of  always  weighing  them 
equally  as  was  done  in  this  study. 

In  an  experiment,  considerable  time  can  be  saved  by  transmitting  a  broad¬ 
band  pulse  instead  of  individual  frequencies.  It  would  thus  be  useful  to  investigate 
whether  this  would  yield  the  same  results  as  a  stepped-frequency  sweep  for  specu¬ 
lar  scattering.  However,  for  extended  scatterers,  the  interference  between  widely 
separated  points  on  the  scatterer  may  not  occur  at  a  given  frequency  (or,  at 
least,  may  not  occur  in  quite  the  same  way)  if  the  wave  does  not  remain  at  that 
frequency  for  at  least  a  certain  length  of  time. 
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6.4.6  Future  Work:  Scatterer  Locations 

In  this  study,  we  have  sought  only  to  find  the  scattering  center  rather  than 
the  location  of  the  individual  point  reflectors  or  the  boundaries  of  an  extended 
scatterer.  It  would  be  useful  to  investigate  how  knowledge  of  the  scattering 
center  would  be  used  in  an  actual  application.  For  example,  could  a  mobile  robot 
avoid  obstacles  knowing  only  the  location  of  their  scattering  center?  If  not,  some 
algorithm  could  perhaps  be  developed  which  could  estimate  the  boundaries  of 
an  extended  object  once  its  identity,  scattering  center,  and  orientation  had  been 
determined. 
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Appendix 

Fourier  Derivation  of  the  Holography  Equation 


A  method  for  the  derivation  of  the  holography  reconstruction  equation 
using  a  Fourier  transform  of  the  wave  equation  was  suggested  by  Eschenberg  and 
Hayek11  and  is  shown  in  greater  detail  here.  As  in  Chapter  2,  the  two-dimensional 
Helmholtz  wave  equation  when  the  pressure  has  a  harmonic  time  dependence  is 
given  by 

V-R{x,  z)  A  k-R{x,  z)  =  0  ( A.l) 


with  the  Laplacian  defined  here  as 


(A.2) 


After  multiplying  all  terms  by  (l/\/2n)exp(—jkxx)  and  then  integrating  over  all 
values  of  X  we  have 


where  J-  and  R  represent  the  forward  Fourier  transform.  We  wish  to  modify  each 
term  of  Eq.  (A. 3)  so  that  it  is  expressed  only  in  terms  of  the  Fourier  transform  of 
R.  This  is  a  trivial  matter  for  the  last  term  since  kx  is  independent  of  the  variable 
of  integration  x.  The  order  in  which  the  integration  and  derivative  are  performed 
in  the  second  term  can  also  be  interchanged  in  a  trivial  manner  since  the  Fourier 
transform  has  no  dependence  on  z.  However,  the  first  term  of  Eq.  (A. 3)  is  more 
complicated  since  the  Fourier  transform  has  a  dependence  on  x. 


A  series  of  expansions  can  provide  the  needed  substitution.  We  begin 


L'.'iS 

by  noting  that  the  second  derivative  of  the  Fourier  transform  with  respect  to  x 
equals  zero  because  there  is  no  dependence  on  x  left  after  the  Fourier  transform 
(it  has  been  replaced  with  a  dependence  on  kx).  This  quantity  can  be  expanded 
as 
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(A. 4) 


or 
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J  kx  R(kx ,  z)  —  0 


(A. 6) 


or 
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Substituting  Eq.  (A. 7)  into  Eq.(A.o)  we  find 
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'd2R(x,z)' 

dx- 


=  -k]  R(kz,  z) 


Finally,  we  substitute  this  result  into  Eq.  (A. 3)  and  use  k~  =  k^  +  kz  to 

d2R(kr,z) 
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hich  is  the  same  as  Eq.  (2.6). 


(A.8) 
arrive  at 
(A. 9) 
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